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ABSTRACT 


This report pi*esents, in two volumes, a recursively formulated, first-order, 
semianalytic artificial satellite theory, based on the generalized method of 
averaging. Volume I comprehensively discusses the theory of the generalized 
method of averaging applied to the artificial satellite problem. Volume 11 (to 
be published in early 1978) presents the explicit development in the nonsingular 
equinoctial elements of the first-order averaged equations of motion. The re- 
cursive algorithms used to evaluate the first-order averaged equations of motion 
are also presented in Volume 11. 

This semianalytic theory is, in principle, valid for a term of arbitrary degree 
in the expansion of the third-body disturbing function (nom*esonant cases only) 
and for a term of arbitrary degree and order in the expansion of the nonspher- 
ical gravitational potential function. This theory has been implemented in the 
Goddard Trajectory Determination System (GTDS) Research and Development 
(R&D) version. 
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SECTION 1 - INTRODUCTION 


In the past, considerable attention was focused on the formulation of the equa- 
tions of motion for complex dynamical problems and on the method of solution 
to insure that a sufficiently accurate result, meeting the investigators require- 
ments, was obtained with an economy of effort. Without such careful consider- 
ation, the most prominent problem of classical mechanics, i.e. , the motion of 
planets about the Sun, would probably not have been solved with anywhere near 
the accuracy actually obtained. It is a testimony to the ability of men such as 
Lagrange, Gauss, Leverrier, Hill, Hansen, and others that not only ingenious 
formulations of the equations of motion were obtained but that the tliousands of 
arithmetic operations required to evaluate the solution wei'e organized in such 
a manner as to minimize the number of these operations and considerably reduce 
the probability of undetected accidental errors. 

The advent of the high-speed electronic conjputer has relaxed this consideration 
by making brute-force, error-free calculations possible. However, the compe- 
tition for computer access has grown rapidly within the last decade. As a result 
of this overload on computer resoui'ces, current problems of interest should be 
formulated in a manner that not only fulfills the investigator's requirements but 
also minimizes computational cost. 

One of the more computationally expensive dynamical problems today is the pre- 
diction and definitive determination of artificial satellite orbits. Maintaining 
reasonably accurate ephemerides for the ever-increasing number of artificial 
satellites (which include active scientific, defense, communication, and weather 
satellites as well as defunct satellites, launch vehicles, and other debris) requires 
a considerable expenditure in terms of computing time. Also, prelaunch mission 
analysis requires that several hundred satellite trajectories over periods of up to 
several years be generated for the purposes of lifetime and geometry constraint 
analysis. In addition, mission feasibility studies corisume an inordinate amount 
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of computer time. Generally, these applications do not require the extremely 
accurate high-precision orbit generation techniques which rely on the expensive 
process of numerically integrating Newton's equations of motion or some equiv- 
alent set of differential equations. 
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1.1 REVIEW OF ORBIT GENERATION TECHNIQUES 

Another approach to the artificial satellite problem is provided by the purely 
analytical methods of solution in which analytical formulas for the coordinates 
or orbital elements are usually obtained to first or second order in a small 
parameter. A standard approach is to separate the short-period, long-period, 
and secular components of the motion through a series of canonical transforma- 
tions (Reference 1). The secular contributions to the motion are evaluated at 
a given time, and the canonical transformation used to remove the long-period 
compcment of motion is inverted to provide the long-period mo^mn in terms of 
the secular elements. Finally, the transformation to remove the short-period 
terms is inverted and evaluated with the secular and long-period contributions 
to the elements, thus obtaining the short-period contributicms to the motion. 

Although computationally efficient analytical satellite theories have been devel- 
oped,^ many of these theories suffer from severely restricted perturbation mod- 
els. Several theories are limited to the lower degree zonal harmonic terms in 
the nonspherical gi\*vitational model of the central body. The third-body pertur- 
bation, when included, is usually restricted to the cases of very cIose-Earth 
satellites. Also, many of these theories are resti'icted further by the ut:e of 
the small eccentricity and/or small inclination approximations. In addition, 
the use of Keplerian elements in these formulations introduces singularities 
caused by vanishing eccentricity and/or inclination. Some of these limitations 
are accounted for by the fact that many of these analytical theories were devel- 
oped manually. The tremendous amount of necessary algebraic manipulation 
required that these theories be severely restricted. 

In ttie last decade, the appearance of machine automated algebraic processors 
has facilitated the development of analytical satellite theories with more sophis- 
ticated perturbation models. All that is required is sufficient computer time 

^Y. Hagihara (Reference 2) gives an extensive list of references to the work in 
artificial satellite theory. 
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and storage. However, a reasonably general first-order .analytical satellite 
Hieory can comprise tens of thousands of terms which require a prohibitive stor- 
age capacity. The <Mily w.ay to reduce the storage requirements for .an explicit 
analytical theory is to restrict the theory itself.^ 

Finally, although several attempts to incorporate atmospheric drag in analytical 
satellite theories have been made, they have proven less than adequate for pro- 
ducing reasonably accurate ephemerides over extended time intervals. This is 
not surprising in view of the fact that even high-precision numerical techniques 
which use sophisticated atmospheric models have difficulty predicting ephemerides 
of strongly drag perturbed satellites over periotls of several weeks (Refercrce Hi. 

The method of averaging offers anotlier appro.ach to the artificial satellite prob- 
lem that has been shown to be more computationally efficient by several orders 
of m.agnitude th.an the high-precision techniques (Itefei'cnce 4). In acklition, the 
method is very flexible with respect to the pcriurbation models and suffers fewer 
I'estt'ictions than purely .analytical satellite theories. Although not .as accur.ate 
as the high-precision techniques, this technique produces results sufficiently 
accurate for all but the highest accuracy requiroments, e.g. , nuaneuvers, etc. 
More specifically, an .application to first order of the method of averaging pro- 
duces the long-periotl and secular motion of a s.atellite extremely accurately in 
moat cases (lleference 4) .and provides for the recovery of probably 90 to 9.*) |ier- 
cent of the short-period motion (Reference 5). Consequently, this apfiroach 
provides a low-cost, long-term orbit prcniiction capability for the following: 

• Mission fe:isibility studios 

• Mission imalysis (lifetime and goonuJtric constraints) 


For certain applications whoiai one particular tjqie of satellite is encountered, 
e.g., circular geosynchronous satellites, a restt'ictod theory is not only accopt- 
:ible but advisable. If, however, a single theory is to be used for several differ- 
ent tj’pes of satellites, a general theory is required. 
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• Trai'king station acquisition schedules 

• Dynamic moileling in iJcfinitive orbit lietcrmination priHvdures where 
eitlter extended data intervals or extemloil dafa gaps are eiu'ounteix'tl 

• Dji’namic motlcling required for differential correction (IH'i pr»- 'dures 
used to solve for dynamical parameters, c.g. , high-order ge * >«c- iai 
coefficients 

The motivation for using the metliod of averagiiig proi'edure is as follows, IV 
maximum step size which can lie used in the numerical integration of a set of 
differential equations is constrained by the highest significant frequency ciwitaim'd 
Aercin. The method of averaging is used to remove high-frequency components 
from the equations of motion- The resulting averaged equations of motion are 
integrated numerically but with a significantly greater step size than can Iw usetl 
with the high-precision equations. The long-jieriod anti st'cular compoiuints of 
the satellite motion are thus obtainetl. The short-|»'riotl compoju'nt of the motion 
can be computet! eitht'r numerically (Mefcrence 5) or from analytical formulas 
which are presimtcti in \ olume II of this report. In most cases, the computatitmal 
savings achieveti by the larger step size (which results in fewer force evaluations) 
far outwieghs the increasei! cost of the derivative evaluation, thereby effecting a 
significant decrease in tlw' overall computational wsts. 

The technique of removing the high-frequency terms from tin' equations of motion 
was first useti liv l.agrange in his investigations of tbt' plaiwtarv motion, ftecause 
of a particular fornuilatiixi of die eijualions of motion tieveloped by I.agrange, the 
high-frequency terms, in the case of conxervative perturbing fori'es, could he 
isolated more or less by inspection. However, a rigorous mxthematical founii- 
ation for this technique vva.- not provided until the relativt'lv ive*mt work by Krylov 
.and Bogoliubov (Ueference 6) on asymptotie methods for nonlinear oseiHationa. 
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Two approaches arc available for the application of the method of averaicing. 
The high-frequency components of the equations of motion can be removed 
numerically by application of a quadrature around an appropriate formulation 
of the high-precision equations of motion. This proc'Cdurc is known as the 
numerical averaging approach. If the perturbing forces are conservative, the 
equations of motion can be expressed using Lagrange's formulation, and the 
averaging quadrature can be performed analytically. Under certain assump- 
tions,^ this method produces the same result as that obtained by inspection. 
This semianalytical procedure of numerically integrating the analytically aver- 
sfed equations of motion is referred to as the analytical averaging approach. 


The assumptions arise when either the Greenwich Hour Angle, i.e., the Karth’s 
rotation, or the fast variable of the disturbing third body apiwar in the pertur- 
bation models. Specifically, these quantifies are assunted to be completely 
independent of the satellite fast variable, both explicitly and implicitly through 
the time. 
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I.? nil': NrVKisu\%t avehmunc. ap»koach 

He^x'ntlv, th*‘ »»vraKiHR niothtxl has Ix'on siu'oodsfully U; tho 

artificial satoUiu* pv«bk'm (Hcfcrcncos -4, a, and 7K This ♦•♦.Unique i« partic- 
ularly flexible in that it can lx? »j,iplivHtto any perturbation »hich can be ^K?tern\- 
inidtically modeled. It is also quite attractite Ixcause of the wase with ♦hich 
d'fferent element sets can be accsifciniotlateii (Hefewnct' Si and Ix'cause of the 
ease of mollifying the first-onk'r assumption (Uoferwnces 5 and^. Numerical 
averaging is also ap{x?aling tweause the implementation of the procedure seems 
to be rather straightfonvard. 



However, the impletm'ntation. »t)d moiv importatulv, the pro|xM- use of numeri- 
cal (as well as analytical) averaging techniques depend on the understanding of 
s*'veral basic concepts, many ofs*hi‘*h are atidresseti in this report and in Uef- 
viv'nces 4 and 0. Furthermore, Early (lleference 10) demonstratwd that a 
straightforward applieation of the numerical averaging technique is not well 
suited loesses where tht' f»fft\irbing force vanes by several orders of magni- 
t«Mk' iner ;♦ short arc id the orbit while remaining essentialli negligible outside 
that interval. 

Nobvithstaiuling, numerical averaging hss lx>en shown (Ketierence l» to Iv an 
ettective procevlure for generating tlw long-periixl xitil sei-ular motion of a satel- 



lite for a w- W variety cases and to Iw i-o«i*iih»raM> more effu-iett tha« the f 

■f 

high-precisiudi techniques. t'lidMeqiwwifly, the numerical averaging appro.U’h ^ 



lii ! 




This is, of course', not peculiar to the averaging method but rather to the form 
of the high-precision equations of motion. 


1.3 Tn£ ANALYTICAL AVEIUGING APPROACH 

The method of analytical averaging is attractive because it is not only signifi- 
cantly more computationally efficient tnan high-precision techaniques but also 
is usually an order of magnitude more efficient than numerical averaging tech- 
niques (Reference 9). This computational advantage is accounted for by the fact 
that the analytically averaged perturbation models, although more complex than 
the high-precision perturbation models, are evaluated only once per integi'ation 
step. The numerical averaging approach requires that the high-precision pertui*' 
bation models be evaluated once at each abscissa of the quadrature. Thus, the 
method of numerical averaging requires between 12 and 96 force evaluations to 
compute the averaged element rates (Reference 10). In addition to the greater 
computational efficiency, the analytical averaging method offers greater pre- 
cision with respect to co 'iputation of the element rates and therefore should be 
used whc"<'ver poss) .e (Reference 8). 

The analytical averaging motliod has been used in the d:;velopment of several 
averaged orbit generator programs (References 11, 12, 13, 14, 15, and 16). 
These programs suffer from one or more limitations, however. In particular, 
most programs are based on theories formulated in terms of the Keplerian ele- 
ments, which produce singularities in the equations of motion fo\’ vanishing 
eccentricities or inclinations.^ Dallas and Khan (Reference 14) modified the 
element set to remove the small eccentricity problem; however, the small in- 
clination problem remains. The Earth Satellite Mission Analysis Program 
(ESMAP) initiated by Cefola (Reference 11) is formulated in a completely non- 
singular element sot but is severely restricted in its perturbation models as 
are the progranis described in Reference 15 and 16. 
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The projiiruni dovoloiwd by Wa^nor (Uoforoiuv i;U is based on nynofnl express- 
ions for the analytleully averajii'd |H>rturbatl(»n nuxlels develo|x'd by Kanin (Uef- 
orenees 17 and I Si whieh are formulated in terms of sinuMlnr Keplerinn elements. 
Cook (Uefereneo 1(11 implemented Kanin's ix'rtnrbation models using Allan's re- 
cursive algorithm for the inclination functions and a recursive algorithm for the 
Hansen ctx)ffii ients based on the recursive propoj ites of U^gendre polynomials. 
Unfortunately, Cook's program 1s baseil on the singular Keplerian element set, 
and the nonspherical gravitational potential is restricted to the zonal harmonics. 

Hxamples of computer-generntoil, explicit analytically averaged |K‘rturbation 
models ace given by Sridharan and Kenard (Heference 101 for the long-|H'riod, 
disturbing tliird-luHiy model using U»e potentially singular Keplerian elements 
and by t'olllns (UeU'rence 201 for a restricted 2:1 resonant gx'opoleni ial nuxlel 
using the nonsingular equinoctiid elements. 
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1.4 UECKNT DKVKLOPMENTS IN ANALYTICAL AVKUACllNG 'lUKOHY 

Very rectntly, several authors have investijjntotl general, analytically averaged 
perturbation models for the third-body and noasphorieal grtvitational perturba- 
tions in terms of nonsingular elomant sets. Cefola agd Broueke (Heferenee SI) 
developed i*eeursivoly formulated models for the nonresonant tliii'd-biHly and /.onfil 
hai'monie perturbatioos based on the equinoctial elements. The development of 
the zonal harmonic model is similar to that of ("ook's model, with the exception 
that the inclination function is developed in terms of associatoil Ixigendre poly- 
nomials and their derivatives and certain complex polynomials, ('efola's thiixl- 
body model is developed in terms of the direction cosines of tht disturbing third- 
body position vector, which proves computationally efficient but is limited to 
nonresontnt cases. Cefola outlined an extension of tiis zonal hnrmonic nuxlel to 
include the nonresonant tesseral harmonic terms (Hefeix'iu'c 22) and later c>om- 
ploted and extended the motiel lo include resonant phinomena (Hflfeixmce 2d). 
Ciiacaglia (Reference 24) rtformulaied Kaula's v>erturbation models (using Allan's 
inclination function) in a nonsingular element set and provided a set of recursive 
algorithms for computational purposes. I'inally, Nacozy and Dallas (Hii'ereiU’e 25) 
also reformulated the Kaula goopotential model (using Allan's inclination fiuu'tion) 
in terms of a nonsingular element set. No recursive algxu'ithms were provided. 

The relatively simple recursive algorithms of ('ook, t'efola, anil Ciacaglia arc 
appealing in view of the alternative of evaluating tlie complicated polynomials 
found in the work of Nacozy mul Dallas. However, (he brute-force implementa- 
tion of recursive algorithms can contribute to computational incfficieiu'y and 
can possibly introduce artificial singularities (not in the equations of motion, 
but in the model evaluation). To insure against this possibility, careful consid- 
eration must be given to tl»e onltring of the terms in the models such tJiat the 
ixH'ursion formub’ proceed in tlie pro|x'r direction to avoid small divisors and 
the amount of recomputation and storage requirements are minimized.* The 

Uvfola has considered tl\e question of tlie efficient iiunjk'iucntation of his theory 
(Heferenee 21). 
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TMa twfKX't iii Ml <Mt«rMrth erf a aai'ias of t«i*i Malfttoiimte afith t)to objootivo of 
4«)|il»muttUng iu th«' Ciodilard Tra^otorr IH'torminatUm Systk'm lUTOK) tU>A«aroh 
aiaf Havi'Kifjmi'ttt (U(bl» rovsitw a .sot of I'aour.sive'ly forwidatod, flr.sf-ovilar ta- 
id)itU'»liy avi'va^Hi oquntiotvi of wottonfor mi arttfioi«l .satollit* ix'vhu'had by 
noarosonaat thhHl-bixi>- muI nouji^ihortoal K»‘.svit«tional la'rturbattons. Thi.s anu~ 
Vftionl .six'rajiuig o.ntinl)itUy ouliaiux'.s tJu' trn).s numorioal avor;iRin>j I'sn^abillty 
(ttefomu'c 4) and provJUos for enitimal .svorjqjy'd pi’rhjrbatioa modoUs for oaoh 
i^aK'lfU' ty|« of (Mt'tvirb.stitHi (Wilii Ibio posaiW*' oxocptlonof tbiril-lnHfy ro8on:ini'« 
oaa«a, wbtoti woiv oottvSiik'rviU. V'M'tial ro.suU.s oht.sini'd for .some trf Hi*' 0|>- 
tiinat avoraiti'd fiorturbation tiiodola tniiTD.s havo Iwon pix’st'iitod in Ri'foranoo >. 

VofofaVs .irorsift'od porturbation nuHtola aiv itdoplt'd for t\w nonro.souant third- 
hiv\v and zonal harmonlo ij»vturb.sHi«s. ’llio (•'.s.sor.at h.nrmonio nioihtl w.sa 
dovolopoil uaiil)!; th*' ainiroM'h outliwd by I'ofola in Rofoi'oiu'o 22. Tho nuxkils 
«tin‘olo(K'd nvro gi'noralizod to handlo I’otroKt'aiV as w*'U a diroot oqutnootial 
%'lonionts tao*' .Viqrftidix .K). 

A.H part o1' tlrta invoatiK'ation« a fairly tWt,ail*'d oouipariaon of tho thooi-ioa of 
I'ofola and (liaoatt'Ua wa.s p*'rfonii*'d. Hrk'fly .statod, tfw tht'orios w»'r*' found 
M hn ba.sioally oqutvalont. Minor difforonooa in t.h*' Hu'orioa InoUulo diffoivut 
nMnaingular oloniont .sota and diffor*'nt oomputational pnvoilutv.s for tha iiwli- 
nntion fum'tion. .Ar^i^uim'nt.s oan Iv ni.%b' ooiu'vrning' tba rolativa aifvantatp'.s 
aMl di.sadvantaK*'.s of Uh'.so nonaiiiii'ular alvniont .srta, but in ro^ard to the ivmovat 
of tlw .siuKularitioa from tho oquation.s of motion, both aiv aoooptablo. itiaoaifha 
oompuh's the ontiiv im'lination funotiiw ri'ouraivoly, roquiring .a moiv oonipli- 
oattd roour.sion ivlation with moro back vahwa of Hi*' funv'tion. I'ofola us*'a 
iWiHiraion fornmlaa for .sororal quantitira oompriaiuK tW imdinatiew futK'tion. 

Ttw iwmraion ia'lation.s aro aimplor, roquiriug fowor bai'k rnluk'a, butmoiw 
t'oour.sion formnlaa aix» mwiWd. 
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Regarding the iinplnmeiitntieB in the GTDS U&D version of the resonant tesseral 
harmonic model, it was fait that this otpabiltt.v should Ite very flexible with re- 
spect to the specific resonant harmonic ti>rms used. TIk' existence of a resonaiu'c 
dictates which terms in ttie potential oapansion are significant to the long-ix'riod 
motion. Know 4ge of the comnion charaeteristiee of these terms and the propei' 
use of the recursive algorithms would have provitled a means for further optimi- 
zation of this model. However, the procedure would lieve been automatic, witli 
the program expecting a certain set of terms. Therefore, fer the purposes of 
flexibility and at some additional computational costs, the eonlidbutions from each 
spherical harmonic term are computed entirely indeixindently from all other terms.* 

Due to the extensive new software for the analytical averaging capability as well 
as to the extensive modifications required to the previously implementeti averag- 
ing softwai'e (particularly the input processor and initiali/alion procedures ;md 
Ae attendant added complexity of executing the GTDS R&l) averaging i-apabilitv), 
it was decided that a system description and user's guide for Ihe (ITllS R&D aver- 
aging capability m'ould be issued under a separate cover. In addition, a document 
extending the numerical results Ix'yond those' presented in Ke'ferenc'c is also in 
preparation. This document will discuss the eomputatioual c()sts in terms of 
machine processing time, the accuracy of tlie analytical avt'raging methods, and 
the proctediu'o and algorithms used to develop an automatic truncation c'apability 
to furtlwr optimize the perturbation uukIcIs for eae'h parlieular ease. 

The current report consists of two volumes. The (l)cen'\ of the method of avi'r- 
aging ia discussed in Volume 1. Volume II presents the ('xplieil ilevelopment of 
a semianalytieal artificial satellite theory based ou (he metlunl of averaging. 

Volume 1 presents a fairly comprehensive diseussi<'n of (lie application of (he 
gonnraUKed method of averaging to (he artifieial satellite problem and (he result- 
ing formulation of thn nvorngt^d equniioas of motimi. in vSei'tion 2 , a discussion 

^The capability to automntietUy select tlie resonant terms was implenu'iited in the 
(5TIKS RgD vorsioa. no speH‘i«I rolntionship lunong fliciu ivS nssuinod. 
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of tlii VtrlattoB of Paramo tors (VOP) formuUtioft of ths equations of motion, 
span which the method of averaging is based. Is presented. Section 3 discusses 
the application of the method of averaging to tiie VOP equations of motion. The 
ovitsrion for the selection of short<-period terras ie discussed in Section 3.1, 
and the geoeralized method of averaging is applied to the VOP equations for the 
caee of a single perturbing function in Section 3.2. A discnssion of the application 
at the method of averaging to the caee of two or more perturbing functions is 
presented in Section 3.3, followed by a description of the modification required 
for tiM s|)pUcation of the raettiod of averaging to cases involving resonance phe- 
nomena in Section 3.4. Next, Section 3.5 addresses the a|iplication of higher 
order averaging theories. Finally, a discussion of the first-order short-period 
variations ia tiie elements and their application to oeculating-to-mean and mean- 
to-osculatiBg element conversions is given in Section 4. 

Volume II presents the mathematical formulation of the aonspherical gravitational 
and nonresonant third-body models required for the first-order averaged equations 
of motion. In this volume, the nonspherical gravitational potential is developed 
in the nonsingular equinoctial element sot, and the zonal harmonic model, the 
combined zonal and nonresonnnt tesseral harmonic model, and the resonant tos- 
aeral harmonic model are isolated. The nonrosonant third-body disturbing func- 
tion is also developed in equinoctial elements and in the direction cosines of the 
third body. All models are pt'osented in what Is considered to l>e an optimal form, 
taking Into account the mlgimisation of the combined computational and storage 
costs while avoiding computational singularities. It is this final form of the 
models Umt was implemented in the GTHS U&I) version. 
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SKCTION 2 - riUO VAHIATK^N OK PAUAMKTKUS (VOP) 1 Ql'ATlONS 

Clnsfiioally, the Variatioa of Parameters (N’l'JP) forimilation of the equatic...s of 
fhotion was usoti to investigate tlie long-period anil secular motion of the planets. 
TIk' VOP formulation was introduced by Kulcr while investigating the mutrial 
perturbations of Jtniiter and Saturn and was later gi'ncrali/ed and completed bv 
Lagraitge (Reference 2GI. Since the primary objective of the curivnt investiga- 
tion is the development of an efficient orbit generation mtthod for the pix'diction 
of ttie long-period and secular motion of artificial satellites, the \'OP fonnula- 
tiott was used. 

In Uds section, a derivation of the basic \’t)P equations is presented in an attempt 
to provide some background information to the reader who is not already familiar 
with the metluxl. Altlunigh the derivation pix'sented is not the most elegant, it 
serves the purpose of explaining tlu' basic principles of the mctluvl and provides 
a logical foundation for the form of the \'OP equations used in this investigation. 
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2. 1 PRINCIPLES OF THE VOP FORMULATION 

The VOP formulation of the equations of motion for a perturbed dynamical 
system requires that the solution for the corresponding unperturbed system be 
known. The unperturbed dynamical system associated with the artificial satel- 
lite problem is the classical two-problem of celestial mechanics. As a starting 
point in the development of the VOP formulation, the differential equation of 
Newton describing the perturbed motion of a satellite relative to the central body 
is considered, i.e., 


k^(m ~ (2-1) 

where "f and r denote the satellite position vector and its magnitude, r is the 
velocity vector, k is the Gaussian constant, m and are the masses of the 
central body and satellite, res|jectively, Q is the perturbing acceleration vector 
caused by conservative and/or nonconservative perturbing forces, and t is the 
time. For mg « m , the satellite mass can be neglected. 

For the unperturbed problem where Q = 0, Equation (2-1) reduces to 

» 0 ( 2 - 2 ) 
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A solution of this system of equations requires six constants of integration. 

These constants are denoted by aj (where i = 1, 2, . . . , 6) or by the vector a* 

The constants are identically the components of the initial position and velocity 
vectors or any set of six independent functions of the initial position and velocity. 
The solution of Equation (2-2) is denoted by the vector function 0o(a, t) . The 
method used to obtain this solution is discussed in References 27 and 26. The 
eolutioB ^0 describes the motion of a point on an ellipse at a particular spatial 
orientation with the central body located at one of the foci. 
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In the VOP formulation, the perturbed two-body problem represented by Equa- 
tion (2-1) is assumed to possess a solution 4 > of the same form as the function 
with the single exception that the constants of the unperturbed motion, a^, 
vary with time. Solving Equations (2-1) then reduces to determining this time 
dependence. 

The VOP equations of motion consist of a set of six first-order differential equa- 
tions .as fol!.,;v8: 


n 0 


\\ i! 

0 


* G,j(a;,t) 


(k - 1,2,...,6) (2-3) 


where the constants of the unperturbed motion, referred to as elements, are 
treated .as time-de|jendent parameters. This system of equations can be obtained 
directly by transformation of Equations (2-1). Eitpressing the three coordinate 
variables in Equations (2-1) formally in terms of the six elements and the time 
results io the throe equations 




(i - 1,2,3) 


\ 


involving six unknowns .aj,. Consequently, throe arbitrary relations or constraint! 
mty bo imposed on the six elements. These relations may be sjecified implicitly 
•nd are usually chosen !uoh that the following equations are satisfied; 


(t‘ 1,2,3) 


_1 _ 1 . 


1 TT 


t < I 




which requires that 


V' Mj ^ 
Z, K dt 


(i — 1) 2( 3) 


The metivetton for this particular choice is discussed below. 

The implicit relations betweea the position and velocity and the six unknowns a|^ 
spectfied by Equations (2-4) end (2-5) will be used to transform Equation (2-1) into 
Equations (2-3). Differentiating Equations (2-5) with respect to the time yields 


* Z-i ba..bi 


bi <Xt 


(i = 1,2,3) 


Substituting the right-hand sides of Equations (2-7), (2-5), and (2-4) into Equa- 
tions (2-1) yields the following three first-order differential equations in the six 
unknowns aj^ ; 


da. 


hi dt 


+ k^m 




“ (2-8) 


(i- 1,2,3; j- 1,2,3) 


Equations (2-6) provide the three other first-order differential equations required 
to determine the system. 

The function fj , representing the ith component of the position vector, is deter- 
mined from the formulas for elliptic (unperturbed) motion, i.e., through 0 q, 
which relate an instantaneous position to a set of instantaneous ulements (in fact, 
infinitely many). It is not imm -diately obvious from Equation (2-4) alone that the 
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perturbed fvlocify vector can ba rslabtd to the aame set of iastantaneoue cletneate 
tturough th«M formulas. However, Equation (2-5) tndlcalea that the velocity com« 
ponents are detemined by dlfferentiatlDg the position functions, f| , while holding 
the elements constant, which la exactly the requirement for unperturbed motion. 
As a reealt, at any time t , the perturbed elements always correspond to a set of 
unperturbed elements. Such elements are referred to as osculating elements. 

The three constraints imposed on the elements by Equations (2-5) are not the only 
set possible, but they are the only set that allow both position and velocity to be 
related to these perturbed elements through the formulas for elliptic motion. 

la Eqaattons (2-3), five elements can be chosen such that they complet«>ly specify 
the osculating ellipse in space. The sixth element, , In conjunction with the 
time t specifies the position of the object on the osculating ellipse at time t . 

The function a, t) represents the time rate of change of the ith osculating ele- 
ment caused by the pertiu'bing force, la most cases, the perturbations are small 
compared with the central force, and, therefore, the magnitude of the function 
Gjj is small. Consequently, in most problems the elements aj^ are sl<mdy var>’- 
ing. 

For conservative pj>rturblng forces, the osculatiag clement rates can be repre- 
sented in terms of the partial derivatives of a disturbing function. The disturbing 
fbnetion is the tegattvo of the potential func'tion, honoe the restriction to (xmeerv- 
ative portujibing forces. To obtain a formulation dcpendeat only ou the elements, 
the disturbing fulction is developed in terms of the olomeots through a formal 

Fourier series expansion. Also, the Fourier series representation petanits 
% 

isolation of specific ft'equsncies in the moMon by inspeotiem. If the series expan- 
sion is developsd literally, Equations (2-3) can be integrated tet'm by tem using 
tte methoil of successive approximations to obtain an analytical approximation 
to the solution (deference 2). This approach ts known as the method of general 
pcrturl)atlons. 


5-r» 


tinier the category of apccial |>erturbaUon tnethoda, aeveral wuncricat tech^^ 
niquea luive been developed for evaluattikg the osculating element rates given by 
Equations (2>3). A particular solution for these equations is then generated 
using a numerical integration procedure. There are essentially two formula* 
tiooo of the special perturbation technique associated wtih the VOP formulation 
of the equations of motion. One formulation, associated with the name of Gauss, 
uses closed form expressions for the osculating element rates, i.e. , the func* 
tiona G}^ are formulated in terms of the components of the acceleration. The 
other formulation is based on a Fourier series expansion for the disturbing 
fuactioii as used in dm general perturbation method except that the coe0icients 
are generated by some numerical scheme. 
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2.2 TJIE GAUSSIAN Vt>ri:<H'AT10Ni 

Thi* particular form of t)w VOP«qu»«i«»* i# amiSfi as^'llows. d'irsf, 

Equation (Z~5) ia aotMtituted into Uquat4»na (2-8) to yi«U1 


a‘f, y- 6ii . , » «i 

(W 

j*i 


•, (i-i.2,;n 


Clenrly, thu corroa(x>tuiiiH: «'quation for the iiftiTTiniMiBninUnn is 


> \»/A 


(i ' l,2,;i) 


8uk)ti'9H'ting Kquation (2-!0) from EqnnhMi (2-9) Kin'O-s 


L da*, • 


(i 1.2.8^ 


Multiplying both aides of b;quations (2-11) by and summing over the index 
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dxi da^ 

•«3l 



(j = 1,2,...,6) (2-13) 
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where is the classical Kronecker delta function since the elements aj^ are 
mutually independent. Consequently, Equation (2-12) takes the form 



(2-14) 


or, more simply, 


L-' 


I ' 

L : 


CL. 



(j — 1* 2, • . « , 6) (2-15) 


This result is known as the Gaussian form of the VOP equations of motion. 

i i 

The right-hand side of these equations can also be formulated in cylindrical 

coordinates where the radial, transverse, and normal components of the accel- i | 

eration are used. This particular form of the equations can be found in most 

celestial mechanics references (e. g. , Reference 29). The Gaussian formulation j | 

is particularly attractive becuase it is appropriate for both conservative and 

nonconservative perturbations. However, because most accelerations are t •«'t ) | 

formulated in terms of position or position and velocity rather than as a Fourier 

series expansion, periodic phenomena cannot be isolated from the acceleration j | 

model by selecting the appropriate terms by inspection. Therefore, a numerical 

procedure must be used for isolating specific frequencies in the motion. j | 
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Because of tte nexibiUty aad relative ease of in-ple^eatatioa. the Gaussian for- 
urulation has been used in the development of numerical first-order averapng 
procedures (References 4, 5, 11, 12, 13, and 141. This foimuiatioi. ins 
disadvantage that conversions from the elements to posibon and voloiMty mu t 
L af^licdlmiver U.e element raUiS are evaluated, i.e. , at every Integra 
Ship, in the bagrangian formulation, this particular disadvanta^ is avoided t 
Jpossible exf^nse of the closed-form expressions for the eguations c motion. 
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2.G THE LAGRANGE PLANETARY EQUATIONS 

The derivation of the Lagrange VOP equations of motion (referred to as the 
Lagrange Planetary Equations) is identical to the Gaussian formulation through 
Equations (2-11), with the exception that the perturbing function or acceleration 
component, Qj , is restricted to depend only on the position and can then be ex- 
pressed as the gradient of the disturbing function, R(Xj, Xg, X3) , i.e. , 




(i = 1,2,3) 


(2-16) 


Equations (2-11) then talce the form 


V A 


(1 - 1 , 2 . 3 ) 


Multiplying Equation (2-17) by 3x./Saj and summing over i yields 


(2-17) 


3 (> 


^ ba.\ “-k ■ 2_, aij 




6R dR 

^ ^ ‘ (j ” 1,2, •••,6) (2-*18) 

OX\ oCXj 


Similarly, multiplying Equation (2-6) by axj/aai and summing over i yields 


3 (> 


W ^ • 

Z. L. 41-i i««.K 


k • 


(2-19) 


i»X k».X 


(It should be noted that Xj has been substituted for fj in Equations (2-6),) 
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St^Viraotimg Iqiiatioft (2-19) from CqaatioM (2-18) yields 


i; 


(j l»2,f«««j6) (2—20) 


where 


f/.. 1 - V/^ ^ 


( 2 - 21 ) 


4s called ttie Lagrange Bracket. 

Although there are a total of 36 Lagrange Brackets required for the complete sot 
of equations specified by Equation (2-20), at most only fifteen must be determined 
because 


[ S » S ] 


(2-22a) 


[S, h] - - [d^, dj] 


(2-22b) 


These conditions follow from inspection of the definition given by Equations (2-21), 
It should be pointed out that the Lagrange Brackets depend only on flie formulas 
for elliptic motion because 


dOLb dO'K 


^ 6 bf; 

do.^ b(Lu tk. 
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The fifteen necesaary La^aage Bracket! required for Equation! (2-20) can be 
evaluated explicitly in terms of the elements and the system of equations inverted 
to yield a^ . An explanation of the evaluation of these quantities is presented in 
Reference 29. 

An alternate derivation of the Lagrange Planetary Equations can be obtained with 
the aid of the folloiwing relation given by Broucke (Reference 30): 

4 

V * / X 

di-, " '2^ (ft-k.O-j) (2-22 

ariiere the quantity (Sj, is the well-known Poisson bracket and is defined in 
Cartesian coordinates by 





da] 

Sir ^ 


(2-24) 


Hie Poisson Brackets also share the properites of the Lagrange Brackets, i.e. , 

® (2-25a) 

= - (O-j. Q-k) (2-25b) 


Equation (2>2S) is immediately verified by direct substitution of the Poisson 
Bracket definiti<m* 
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Ei{a«MiQg the Gaussian VOP equations (Equations (2*-).S)i in tsims ai the dis- 
turbing function yields 





6R 

da-, 


'-'muiNAL 


OFi?(X)R^QUAlIlTY * * 


(2-26> 


Substituting the eipression for da|^/dX| in Equation <2-23^ into Equation (2-26) 
immediately yields 





do.j 


(k = l,2 6) (2-27) 


or simply 

‘k * 

Equations (2-28) are the Poisson Braclcet representation of the Lagraixge Plane- 
tary Equations. 

The relationship h. '.;;-'en the Lagrange and Poisson brackets is iminediatelv 
obtained by substituting Equation (2-26) into Equation (2-20). The result is 


■I 


/ V 


(k - 1,2, ...,6) (2-28) 


4 4 






(2-29) 
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vhlcii requires the conditUm 


Zj ^■K ] 


(2-30> 


(Equation (2-25) was used to remove the negative sign in Equation (2-30). ) 

The particular VOP formulation adopted for this report is a modified version of 
liSgrange's Planetary Equations and is given by 


\ * , \ 

^ (i = l,2,....5) (2-31a) 


« . 'Tr. 


(2-31b) 


where n is ^ mean motion and Sg now denotes the variable X under the summa- 
tion. The variable X , referred to variously as the fast variable or the rapidly 
rotating phase, is not a true slowly varying element but is a linear combination 
of the time with an element such that 


X ■ fit ^ Olg 


( 2 - 32 ) 


The parameter X meaaures the angular distance of the satellite from aome 
departure point in die orbit. This modification, which was made by Tisserand 
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(Rflf«roQce 31), is nec«s8arv to avoid the presence of mixed seeular terms in 
tim equations ctf motion. A mixed secular term has the form 

t*' Got mt 
ain «ni 

•ad quickly degi'ades the solution as time t increases. The appearance of such 
terms is not inherent to Ae problem but to the formulation of the problem. The 
mean motion, n , enters into Equations (2-31) through Equation (2-32). Use of 
the variable A aiipears to have significantly changed the form of the Lagrange 
Planetary Equations. Hou’ever, the original form of tlie equations given by 
Equations (2-28) is easily recovered by modifying the disturbing function with 
the addition of die negative of the total energ>’ to the original disturbing function, 
i.e., if the semimajor axis is denoted by a, then 


R' - R ♦ ^ 


Equations (2-31) can then be expivssed as 




da] 


(2-33) 


where a^ ia understood to represent Uie variable jI . A more complete discussion 
of this question is presented by Plummer (Uefei'unce 32). This I'efinement ia not 
naceaeary for the put'poae of thia investigation and, accordingly', will not be used. 
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2.4 DISCUSSION OF ORBITAL ELEMENT SETS ^ 

The preceding discussion of the VOP equations has made numerous references ' ^ 

to ttje "elements" or "osculating elonients." The question of which element set i ; 

I ; 

to use has not been addressed, and, in fact, a general discussion of the VOP 
formulation need not be conoemed with any specific element set. However, the i 

application of the VOP equations does require the selection of a set of elements. 

There are several well-known element sets, the best known of which is the set i 

of cl. sical or Keplerian elements. The VOP equations formulated in Keplerian 
elements contain the eccentricity, e , and the sine of the inclination as divisors 
and therefore are singular for vanishing eccentricity and/or inclination. There 
are several nonsingular element sets available, and the choice of a particular 
set is arbitrary insofar as removing the singularities from the equations of 

! 

motion. However, some of these sets can present a slight computational advan- j 

tags over otfier sets when converting from elements to position and velocity. 

1 

For other applications, such as differential correction and error analysis pio- i 

oo<hires, the choice oi the element set may no longer be quite so artibrary. 

According to Broucke and Cefola (Reference 33), the nonsiogular set of elements, j 

which are called equinoctial elements, can possess marked computational advan- 
tages over other nonsingular element sets. ; 

The equinoctial element set, a = (a, h, k, p, q, X) , is used in this investigation. . ) 

It is defined in terms of the Keplerian elements by the following: 

a « 0. .J 

* c am . | 

k * € con TA) 

; ) 

^ coaH. 

A » X in 
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wh0i*4 I is the retrograde factor which takes on the values 

I = 1 (for 0^l4(n/2)) 

I = -1 (for (r/2)<l<-rr) 

k more complete discussion of this element set# includiag the Lagrange and 
Poisson brackets and the conversion to position and velociiy, is presented in 
Appendix A. 

The VOP equations expressed in equinoctial elements take the form 
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The diatnrbing functione presented in Volume II of this report arc better expressed 
in terns ai fce direction cosines ( a , ^ , 7 ) with respect to the equinoctial ref- 
erence frame ("f, g, w) of either the equatorial z' axis or the third-body position 
vector) rather than in terms of the equinoctial elements p and q . Consequently, 
expressions of the form 


dp 


dR dot 
det dp 

dR dot 
dot d(^ 


dR ^ 
djB dp 

dR dft 

dft d«j 


dT dp 

dR dY 
dt d({ 


will be used to modify Equattons (2-34) in order to accommodate the particular 
form of the disturbing functions. The following results, presentod here without 
proof, are demonstrated in Volume II: 


dot a , V 

df • - 1 ♦ f) 


(2-35a) 
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whcr* C ia defined as before, Substituting these eaptvesioos into Equations 
jrtelds the final form of the VOP equations cf motion used in the current investiga-' 
tion^ i« e« • 
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Wkere A* and C ana Bafined as Hk S^aiMion (2-34) and 




for any two variables x and y . 

It should be pointed out that a considerable siwiplificatton ocsurs for tlyp Apii' 
resonant ttiird-body and zonal harmonic pertusBitions whets 



O' 

nUs simplification was first reported in C^natisio (5-57) of Referen#ai and 
will he demonstrated in Vtdume II of this rapaot. 





I 


4 



•14 ^ I ! 

i I I ' ^ 



•1 



.- r ■ .^ 

i ■ 



r ; ^ ■ : 


SECTION 3 - THE AVERAGED VOP EQUATIONS OF MOTION 


Classically, in the investigation of the long-period and secular motion of the 
planets, the Lagrangian Variation of Parameters (VOP) equations (known as the 
Lagrangian Planetary Equations) were expanded in a literal Fourier series and, 
with the proper assumptions, the terms which contribute to the long-period and 
secular motion could be easily isolated by inspection. This technique produces 
excellent results when the perturbations are small, and it has been used exten- 
sively to investigate the planetary motions over long time intervals. 

Alternatively, the long-period and secular contributions to the motion can be 
systematically isolated by applying the method of averaging to the VOP equations 
of motion to eliminate the short-period contributions. The solution of the result- 
•ng system of averaged equations is a set of parameters, usually referred to as 
mean^ elements, that describe the long-period and secular deviations of the 
pei'turbed dynamical system from the unperturbed system. 

The technique of eliminating the short-period terms from the equations of motion 
was without a mathematically rigorous foundation until the relatively recent work 
of Krylov and Eogoliubov (Reference 6) on asymptotic methods for nonlinear 
oscillations. The theory of the method of averaging is based on Poincare's 
theory of asymptotic expansions (Reference 34) and the introduction by Krylov 
and Eogoliubov of the concept of a near-identity transformation. The theory has 
been extended most notably by Mitropolsky (Reference 35). 

Further elaboration and discussion of the theory has been contributed by several 
authors. Kruskal (Reference 36) remarked on the possibility of a recursively 
formulated general inversion of the near-identity transformation. Stern 


The mean elements are defined operationally as the solution to the averaged 
equations of motion. Consequently, the exact definition of a particular set of 
mean elements depends on the interval over which the equations of motion are 
averaged. This is discussed in more depth in Section 3. 1.5. 
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(Hoforence 37) dovelojied this recursivo algoritlu'n explicitly. ' Kyner (Refer- 
ence 38) and J. A. Morrison (Reference 39) have shown the Von '/eipel trmis- 
formation method to be a special case of the (;euerali/ed method of averajjing, 
at least to second order, thus establishing a direct link to the mctliods used in 
developing analytical satellite theories." F. Morrison (Reference -10) has 
presented a lucid discussion of the first-order application of the metliod. A 
discussion of the generalized method of averaging is also given by Nayfeh (Ref- 
erence 41). 

Although the discussion in this section is equally valid for man\ ^ '.hcr dynamical 
systems, the primary objective of tliis report is the application of the methoil of 
averaging to the equations of motion for an artificial satellite. Consc(juently, 
the concepts of short and long period are developed in this context. Also, since 
the method of averaging can be applied to either the Gaussian (Kquation (2-15)) 
or Lagrangian (Equation (2-31)) foimuilation of the VOP equations, Uie giMieral 
expression 


da; 

dt “ 1^2 5) (3- la) 


di. 

dt 


£ F ^( o ,,0 


(3- lb) 


(where a consists of the five elements ) is used in the following discussion. 


This recursive algorithm is a general expression relating tlie jth-order term 
in the noar-idc tity transformation to various combinations of the lower order 
terms in the transformation witli lower order contributions to (he mean element 
rates. This I'ocursive algorithm is quite distinct from tlie recursively formu- 
lated first-order theox'y presented in Volume 11 of this report. 

An analytical satellite theory esm be developed using successive applit'alions of 
the method of averaging to remove first tlxe short-period terms and tlien the 
long-perlotl terms. 
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In this section, the generalized method of averaging is applied to the VOP equa- 
tions of motion to obtain systematically the equations for the long-period and 
secular motion. A discussion of the criteria for the selection of short-period 
terms is presented in Section 3.1, and the averaged equations of motion for a 
s' jle perturbing function are derived in Section 3.2. Section 3.3 extends the 
application of the method of averaging to cases with multiple perturbing func- 
tions. Next, in Section 3.4, the modifications required to extend the application 
of the method of averaging in the case of resonance phenomena is presented. 
Finally, the application of higher order averaging theories is discussed in Sec- 
tion 3.5. 
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3. 1 CRITERIA FOR SELECTING SHORT-PERIOD TERMS 

The criteria for distinguiirtiing short-period terms are, in general, subjective. 
The shortest period of significance in the equations of motion effectively con- 
strains the integration step size. For efficient computation, it is desirable to 
maximize this step size while i*etaining the essential character of the motion 
over an extended interval of time. This is the primary consideration in the 
selection of appropriate criteria for distinguishing short-period phenomena. 

To illustrate this point, the following simple differential equation is considered: 

a. = Cj COS [ j (£ -lo)] 

In general, the minimum number of function evaluations required to integrate 
a function of this type over one period is four. The cosine function has three 
zeroes in the interval of one period. In view of the Fundamental Theorem of 
Algebra, any approximating polynomial which is valid over one complete period 
must be of at least third degree. Consequently, the function must be evaluated 
at four points to determine the coefficients of this third-degree approximating 
polynomial, or, equivalently, the function and its first three derivatives can be 
evaluated at a single point, requiring four function evaluations. This does not 
mean that four function evaluations per period provide the best representation 
of the element rate in the example, but only that this is the minimum number of 
function evaluation*^ per period required to obtain the gross behavior of the real 
solution. The accurate integration of such a periodic function using arbitrai'y, 
equally spaced abscissae would probably require six, and more likely, eight 
function evaluations per period, requiring a corresponding number of integra- 
tion steps. 

A useful criterion for the selection of long-period terms is provided by careful 
examination of the frequencies in the artificial satellite problem. 
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y.l.l Satollite-lX'peadent Freuwm’iwM 

The perhirbini; functions Fj(a, Jl) in Kquations nre assiuuod to be Zn |ieri- 
odie in thB satellite fast viriabte, t • Sonw of the slow variables are angular 
<]uaatiti«6 (kepterian ulemunts) or functions of angular ((unntities (equinoctial 
eleuinnts) that produce fundamental periods in the ntoUon of order 0(e~^) . if 
Pj denotes the fundamental period produced by one of the slowly varying angles 
and if the fundamental {lerioti produoed by the fast variable X is 2;r, tlien the 
fundwnentnl periotl, Pl, satisfies U»o relatitw 

i -TT-i 


If the quaatitv cl F. I « I . the period PI must be such that VV»27r , i.e., 
it is much greater than the periods contributed by terms containing the fast vari- 
able Jt . In addition, the \'OP formulation implicitly nssumes that the quantity 

e IFj I is not large. This discussion suggests that terms ileiwmdent on Hie 
' » 'max 

satellite fast variable X and all multiples of X (i.o., mX • where m l,2,d,...), 
which are of periotl 27T m , be considcretl to be short iieritxlic as cooiparetl with 
terms containing the slowly varying angulai' quantities, ('oiutequently, all terms 
with iH'riods of the same order of magniftule as the satellite period ami all smaller 
periotls will lx> ctmsiileretl to l>e short ix'riod terms. 

Other variables which can introiluce short-iK'rioil effects also apix>ar in the |x'r- 
turbing function. More six'cifically, the effects on the satellite motion eaused 
by tlie fast variable of tl»e tlisturbing third body (i.e., Moon, Sun, etc.) or the 
Oreenwich Hour Angle in the nonspherieal gravitational pott-'ntial model must be 
eons Ide red. 

2.1.2 Thinl-lUxlv Kffects on the Motion 

The presenee of the disturbing third-body fast variable in tlie ecpiations of motion 
will eontribute terms witli a fundamenttd period of approximately 2S clays for the 
Moon and 1 year for the Sun. .'lither of these can certainly be I'onsidcrcd to 


2 - 
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produce long-period effects (relative to the satellite period) in the motion of the 
vast majority of artificial Earth satellites. An infinity of multiples of the third- 
body fast variable also appear in the third-body model. Such terms will contribute 
flie periodicities P* to the motion of the satellite where 


* n 


(n — 1^2^3... 


and where P* is the fundamental period produced by the fast variable of the dis- 
turbing body. 

Clearly, as n increases, the periods P* decrease; therefore, vei'y high har- 
monics in the third-body perturbation model will contribute terms with periods 
similar to that of the Earth satellite, thus introducing third-body-dependent 
shcrt-period terms. However, in the absence of resonance, the coefficients of 
these high-harmonic terms are very small in magnitude, rendering the contri- 
butions of these terms insignificant.^ Consequently, the third-body motion (in 
absence of resonance) contributes significant effects with periods of P*/n , 
where n usually remains a small integer. Such periods are, in most cases, 
still very long compared with the periods of most Earth satellites. 

However, certain classes of satellites (e.g. , interplanetary Monitoring Platform 
(IMP) satellites) have orbital periods comparable to the periods of the lower har- 
monic lunar terms cited above. For this class of satellites, the lunar effects on 
the motion cannct be considered to be long period. However, in the ease of a 
strong resonance, a long-period component of the motion is introduced. The 
period of the resonant or critical term is significantly greater than the period 
of the satellite. 


In resonance, the commensurability between the mean motion of the satellite 
and the mean motion of the third body or the Earth's rotation rate causes the 
appearance of a small divisor in the coefficient of the critical term, resulting 
in a significantly increased magnitude for the coefficient and a corresponding 
increase in the contribution of die term. 
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8.1.3 Nottaphtrlctl Gravitatiopal Kffaots on tfaB Motion 

In tha case of psrturbiog effects cauaad by th« noosphericitv of th« Earth's t;rnvi- 
Utional field, the roUtlon of th« Earth, rttprosantad as the Clroaua'ioh Hour Antt'Ie, 
# , contributes terms with a fundamental period of 24 hours. These terms can be 
conaideiwd long period only for close-£ai'th satellites with periods of at most a 
few hours. These "long-period" contributions of 24 hours and fractional multi- 
ples tttei'eof should nut be grouped witli the long-period contributions of several 
days or more caused by the third body ami the slowly vai’ylng elements of tho 
satellite. Conseipwntly, these nonsphericnl gi'avitational contributions will Ix) 
referred to as "medium-perioti" contributions. 

Multiples of tho Greenwich Hour Angle, m(l {i\\ 1,2, ...1, aiH^ear in the spher- 

ical harmonic eapansion representing tho Earth's gravitational model. The har- 
monics of moderately high degree (i.e., m 11, 12, 13, ote.l will contribute 
tei'ms with periods of s few hours or less. Even for olose-Karth satellites, these 
terms obviously cannot Ix) eousidcrod to he medium |wriod and will be referreil 
to as fl-dofiendent short-period terms. As in tlx) Uiird-bocfv ease, tlie eoeffieients 
of those high-degroo harmonic terms are small, escept in the case of I'esonance, 
and produce little offoet on tlie motioti. Consequently, to first order the harmonii’s 
of lower ilogree can Ihj eonsiik'red te prtxluce relatively insignificant medium- 
period effects on elose-Earth satellites, except in the ease of resonance. For 
satellites with largi'r orbital i»'riods, even the medium -jH'riod effects produced 
by die low-order tesseral harmonics must Ix' consick>rod as short perioii. 

In summary, Uie key to the designation of short-|w'riod anil long-|x'riod tei'ms is, 
of eourae, tlie orbital iwriod of the satellite. All twrioils inlnxlueed by the satel- 
lite f.aet variable are eonsidei'Od to be short periixl and ai'e referred to as aaiellUe- 
(tetumdent or Ji -de(wndent slKM't-pertixl terms. 'IIh' otliev froqix'neies in tlie 
<H-nan\ieal system, i.e., tlie freqiwnciea introduced by tlu.' third body and by the 
rotatioa of tlx) central bixiy, imiat be considered in relation to the natural frequency 
of the satellite. 
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3.1.4 


)ilcationa for the Appllcatloo of the Method of Averaglt 


The method of averaging is best suited to cases whei'e quite distinct groups or 
families of frequencies are present. Each of these distinct families is intro- 
duced by its own source, and the distinction is found in tlie specific frequencies 
and amplitudes introduced. Occasionally, the higher frequencies in one family 
approach the primai'y frequency in another family and the separate contributions 
become more difficult to distinguish. Furthermore, elimination of one of these 
families of frequencies by a single application of the method of averaging does 
not eliminate the similar frequencies contributed by the other family.^ 

Additional applications of the averaging pi'ocedure are expensive in the numeri- 
cal averaging approach or require multiple forms of tlie analytically averaged 
equations of motion necessary for all cases that might be encountered. Also, 
multiple applications of the averaging procedure are not always suitable as a 
technique for developing a reasonably accurate orbit generator. In contrast to 
a second averaging procedure, other means sometimes exist for eliminating 
unwanted high frequencies in the motion. 

Proper restriction of the tesseral harmonic terms in the nonspherical gi'avita- 
tional model will eliminate the t?-dependent short-i>eriod terms tliey introduce 
into the motion. Such a restriction has no effect on the secular motion, at least 
to first order, since the tesseral harmonics produce no secular contributions to 
the motion to first order (Ueforence 2). In fact, for all nonresonant satellites, 
it ia recommended that tlie contribution of all tesseral harmonic terms be de- 
leted from the averaged equations of motion. 


In the case of exact resonance, tivo of the families of frequencies are no longer 
distinct. The frequencies in one of the families appear to be integral multiples 
of the fi’equencies in the other family. Furthermore, a single application of 
the averaging procedure will remove all frequencies contributed by both soui'ces 
up to a cut-off frequency specified in tlie averaging operation. This is discussed 
in more detail in Section 3.4. 
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Tl») inc'lusion of these medium-period and tJ-tle|ien(ient short-period eonti’ibu- 
tions in the evaluation of the mean element rates severely restriids tlie step 
size in the numerical integration procedure. The medium -jxjriod conti'ibutions 
have periods of 24 hours or less and they necessarily restrict the integration 
step size to at most 3 to 4 hours. Altliough the amplitiules of these terms are 
not negligible, they do not significantly affect the long-term motion as c'ompared 
m’ith the long-period and secular contributions of the zonal harmonic's. Further- 
more, these medium-period tesseral harmonic contributions can be evaluated 
analytically in the same manner, ;md at the same time if desired, as tlie short- 
period element variations discussed in Section 4. 

If the medium-period effects contributed by the low-order zonal Icarnionics 
are retained in the equations of motion, tlie tl-de|x'ndent short-period terms 
should still be eliminated as described above, since it is inconsistent to elimi- 
nate the satellite-dependent short-period terms while retaining the -dependent 
terms with similar periods. This, in effect, defeats the whole purpose in the 
application of the metliod of averaging by imposing small step sizes in the numer- 
ical integi'ation proc-edure. The arbitrary imposition of larger step sizes in this 
case causes these (1-dependent short-jx’riod terms to introihu'e spurious noise 
in tlie mean element rates and, c'onsequently, in the numeric'ally integrated solu- 
tion. This is explained by the fac't that the contributions of these short-period 
terms are propagated through the numerical integration as though it were part 
of the contribution of a term with a period approximately six to eight times the 
step-size interval. 

The effects caused by Uie tliird bexly can be c'onsidered as exclusively long-period 
for the vast majority of artificial satellites. However, for very-long-ix'riod 
satellites (such as the IMP class with periods of soveral tlaysi, the tliini-body 
(lumir) contribution can in no way be considered to be long (x'riod and tlie appli- 
cation of tlie methotl of averaging must be reevaluatoil in this light. 
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The usual procedure ia these cases has been to use Gauss' method of secular 
perturbations (Uefei’ence 32), which is also referred to as the method of double 
averaging. In this approach, the method of averaging is applied twice in suc- 
cession, once to remove the satellite-dependent frequencies and then again to 
remove the third-body-dependent frequencies. While this method does isolate 
the secular motion of the satellite quite well, the periodic variations contributed 
by the third body to the motion of the satellite may have amplitudes of several 
thousands of kilometers. The elimination of such contributions is usually not 
suitable for generating a reasonably accurate satellite ephemci j... I hc alter- 
native of using a high-precision technique to generate a satellite ephemeris 
should be strongly considered in this case, since such large step sizes ai'e 
appropriate even for the frequencies in the high-precision case. 

A strong resonance in the problem introduces a long-period contribution to the 
satellite motion of considerably lai’ger period than either the satellite or lunar 
periods. In this instance, a single application of the method of averaging will 
isolate these contributions to the motion. However, due to the strong short- 
period variations in the problem contributed by the fast variables of the satellite 
and third body, a second or higher order averaging theory is probably required. 
This is also probably true for the double averaging approach discussed abo.'o. 

Based on the above discussion, a single application of the method of averaging 
is used in the development of the semianalytical theory presented in this report. 
The ^-dependent short-period terms will be eliminated by appropriate restric- 
tion of the potential model. Although it is not recommended,^ the theory for the 
medium-period contributions to the equations of motion will be develo|xjd. llie 
third-body theory developed in this report is restricted to nonresonant cases only 
and to satellites with periods significantly shorter than tlie third-body orbital 
period. 


The analytical formulation of the medium-period contributions has not been 
implemented in the GTDS U&D version. 
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3.1.6 Mean bllemcnts 



Sinee the mean dements are defined o|MM’ationally as the solution of ;he aver- 
aged eqtlfttions of motion, the exaet definition of a speeifie set of mean elements 
depends on the assumptions or eonstraints imposed in the development of tlw 
averaged equations of motion and on the interval over whieh the equations of 
motion are averaged. In this report, the averaged equations of moiion are de- 
veloped so that the moan elements obtained are, in principle, eo. ivalent to the 
mean or average over the averaging interval of the oscuhdmg elements. This 
is demonstrated in Section 3.2.2. The averaging interval (in the absence of 
resonance) is selected to be the satellite iK'rioil to ensure the elimination of all 
satellite-dependent short-iieriod terms. 

The dependence of the definition of a particular set of mean elements on the 
averaging interval has contributed to some confusion in the communication of 
I'csults obtained by different investigators. For many, the term "mean elements" 
is immediately associated witli tlie double-primed elements obtained by llrouwer 
(Reference 42). This element set reflects only the secular motion of the artifi- 
cial satellite. The single-primed element set obtained by Brouwer in the same 
referenc'e reflects both the long-iK'riod and secular motion of the satellite. The 
single-primed element set was obtained by the application of an averaging ojx'r- 
ation over an interval eciual to the jK'riod of the satellite and, conseiiueiitly, is 
the analog to tlw me:m elements used in this rejHJrt. 

In an attempt to eliminate tlie confusion c'aused by terminology, several other 
names, including single averaged elements jmd long-period elements, have been 
suggested. However, ticese terms do not adeejuately define the elements. This 
is because the mean elements are defined wlcolly by the theory from whic'h they 
are obtained and, thei'efore, no simple naming device can adeciuately describe 
them. To eliminate confusion when comparing separately obtained I'esults, Uk' 
con*esponding theories must be understood. Therefore, the terminology’ "moan 
elements’' will be used in this report, recognizing tlie inherent ambiguity in the 
phrase and also recognizing tlio lack of a satisfactory alteniative. 
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J.2 THF, AVERAGBD EQUATIONS OF MOTION FOR A SINGLE PEKTURBINC 
FUNCTION 

The followiag set of differential equations is consideMd: 


dOL'i 

IT * 


(i = 1,2,. ...5) (8~2a) 


n(aO 


(3-2b) 


whdfe the vector a consists of the five slowly varying elements a. . The nsnr- 
ideMlty transformation from (a, jt) to (a, X) is assumed to take the form 


'* * ”*• •»* (i - 1,2,..,,5) <3-3a) 


^ f 

t • I* 5] t' +ou'*S 


l3-3b) 


where the functions j nre 23T periodic in i . The barred variables are ref- 
eri'Od to as moan elements. The quantity c is assumetl to be a small paraDtstev« 
e.g,, a coefficient in oie of ths terms of tf\e spherical harmonic •xpOBHion oC 
the geopotential model or the ratio of the semimajor axes of the satellite and 
ftird-body orbits in the aeries SKpansion of the third-body disturbing function. 

The presQAce of such a small parameter is basic to ths method of averaging. 
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ll the appllcatka of the metliod of av«r»giag, the transforin of the origiiuil ejreleib 
of equation* (Equatioa <3*2>t 0«<«t the oquatioiui of moUuii for the meait elemeotst 
is assumed to be of tiie form 


dil\ A . 

^ ♦ 0(t ) (i * • • • 5) (3~4a) 


4I 


it * ♦ OU***! 

i«i 


(3-4b) 


so that the rate of change of the mean elements clepeikls only on the alowiy vary- 
ing mean elements. 

Basically, the procedure for obtaining the mean elenient equations of motion ie 
to express both sides of Equation (3-2) in terms of the mean elements (f, 4). 
Equations (3-3) and (3—1) are used to transform the left-hand side of Equation (3-2). 
The pt;rturbing function on the right-hand side is expanded in a Taylor series 
about the inean elements and then rearraugeil as a i^ower series in the small pa- 
rnnieter e . The resulting equations are a\o raged such that all tk|a?ndent^ on 
the moan fast variable is «limin<ite<l. The final result yield* onler-by-order 
o}iprossions for the mean element rates, |( a ) » in terms of auitably nveraged 
functions of the perturbing funetion and its partial derivntlven. 
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3.2.1 yorWttUticm of the VOP Equations in Mean Elements 

(,V2» are cxpre<!t«<l in terms of mean elements as folk»i's. Vi»sV 
Eduation (3-3) iv< differentiated, obtaining expressions for the OK^tating efc»- 
menl rate* which depemi only on the eorrespondiag mean elements and fheee 
ratoe. i>«.. 
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0 1,2,... .5) (3-5a) 


di 

dt 






(3-5b) 
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wlwre is »mderstood to designate 1 umler the summation. Substi'utia* 
Edualions (3-4) into b’nuations (3-5), thus introtfiK'ing tl>e fuuotions into 

fV « 4 |tt;tiions of motion for tht‘ osculating eU'ments, results in the expression* 
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The form of these expressions can be rearranged to give 
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+ n(10 


(i = 1,2,...,5) (3-7a) 
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k>l p<l 


n(aO + Afe.i Ct) + yilaO 
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k«l P‘k 


(3-7b) 


The summation over p is not performed for j = 1 and thus does not contribute to 
the first-order terms. 
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Next, the perturbing functions on the right-hand side of Equations (3-2) are 
expanded via a Taylor series about the mean elements as follows: 
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(i - 1,2,...,6) (3-8) 


Q. s S. 
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where Aaj^ - aj^-aj^ are defined by Equation (3-3). The notation S/(Saj^) 
denotes the operation 

^ m 5.^ 

and for the sake of conciseness will be used throughout this report. Rearranging 
Equation (3-8) as a power series in e yields 




(i - 1,2,...,6) (3-9) 


where 
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etc. The moan motion, n(aj) , is also expanded in a Taylor series al'out the 
mean clement, ifj , i.e. , 
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Uearranging' Equation as a power vseries in < yields 


«(a0 *y t*'N^Ca,l) 0(4***M 




where 




(.-i-ldai 


'-if, ’Zu 

nA'O- 






Njd.l) * * it fT ’Zi,! ^ e I 


:i-i7 


I 


1 

( 


] I 


i 





etc. Substituting Equations (3-7), (3-9), and (3-12) into Equation (3-2) completes 
the transformation. Equating terms with like powers of e yields the expressions 
for the jth-order contribution to the osculating element rates, i.e. , 
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(i 1,2,.. .,5) 
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3.2.2 Eliminav io a of the Fast Variable Dependence 

In order to determine the averaged equations of motion (Equation (3-4)), the func- 
tions A. . , which depend only on the slowly varying mean elements, must be 
J 

related to the perturbing function or its power series representation. At first 

glance, it appears that this is accomplished in Equations (3-14). However, the 

functions T) . . are as yet undetermined, except that they are constrained to be 
M J 

2tt periodic in the mean fast variable, I . Fortunately, this condition permits the 
elimination of the mean fast variable dependence. Integrating both sides of Equa- 
tion (3-14) over the mean fast variable, I . on the interval [o, 27 t] eliminates 
the function arj /Si . This procedure of definite integration is referred to as 
the averaging operation and is written as 
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Some properties of the averagUlg ojieration derived from the above definition ai'e 
as follows. If X( a, JL ) and Y( a, i ) are two functions (appropriately continuoua 
and differentiable) which are 2tt periodic in JL , then 
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X(S,I) Y(S.,i) > I* 
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X(a,I) Y(l,ly » (3-i6c) 


/oX(oL,i)^ = yo ^X(l,i) 
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(3-16d) 


(k 1,2, ...,6) (3-16e) 


where p is any function independent of X. These properties will lx) used implic- 
itly throughout the remainder of Uds section. Ik'cause 10. . is 2n periodic in £ 


(a condition of tlie near-identity trjmsformation). 
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In view of Equation (3-17), the averaging operation also yieUis 
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As a result, the averaged equations representing the jth-order contribution to 


the mean element rates are 
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These equations can be simplified by requiring that 
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which follows from the application of the averaging operation to Equations (3-5). 
Consequently, the mean elements (a, I ) represent the long-period and secular 
contributions to the osculating elements (sT, i ) to within a constant, and 



(3-22) 


where • is a constant. Equation (3-22) follows from Equations (3-20) and 
(3-17), A logical extension of the constraint in Equation (3-20) is to require 
that these constants vanish identically, i.e.. 
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Initially, in the development of the averaged equations, the functions Yl. . were 

M J 

quite arbitrary except for tbc condition of 2 periodicity in i . Equation (3-20) 
restricts these functions to contain only short-period, mixed short-period,^ and 
constant terms. Equation (3-23) further restricts these functions to pui’e and 
mixed short-period terms only. That such restricted functions can be deter- 
mined is demonstrated below. 

Applying the constraint expressed in Equation (3-20), Equations (3-19) reduce to 



(i = 1,2,...,5) (3-25a) 


A mixed short-period tcmi is the product of a piu’e short-period term, g(I), 
and a long-period term, f(i) , i.e. , f(S) g(i) . 
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The averaged equations of motion are now completely specified in terms of the 
expansion of the perturbing function and, in the case of the variable i , the ex- 
pansion of the osculating mean motion. More explicitly, substituting Equations 
(3-25) into Equations (3-4) yields the following expressions for the averaged 
equations of motion: 


0 * * 1 ^** 




(i - 1,2, ...,5) (3-26a) 
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(3-26b) 


The functions fj and for k > 1 ai'e formulated in terms of the as yet 
undetermined short-periodic fiwctions ^ . This dependence is sho^»n explic- 
itly in Equations (3-10) and (3-13). The averaging operation does not fi'ee the 
averaged equations of all contributioim from the short-i«r iodic terms. Such 
contributions are, in fac*t, the source of the higlier order terms in the averaged 
equations of motion. The product of two short-period functions can yield a long- 
period term; for example, in the product 

Jh(i) »ir»i ] nini] • q(a,) - ^ 

both factors are of short period, yet the protluct contains a long-period term 
(i.e., a term independent of i ). 
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bispectioa of Fquatioos (3-10) and (3-13) indicates that products of the partial 

derivatives of the osculating force function, Fj , with the functions tl- . appear, 

‘ *» J 

as do products involving two or more of the functions fli, j • Such products can 
produce long-period terms, as described in the above example, that will always 
be of second or higher order in the sm ill parameter. 

3.2.3 Determination of the Short- Period Functions, T^j.j 

The general formulation of the averaged equations of motion is completed by ob- 
taining the functions j from the information contained in the method of aver- 
aging. In the following discussion, these functions are determi'»''d without the 
constraints expressed by Equations (3-20) and (3-23). How'ever, the justification 
for these constraints is demonstrated. 


A partial differential equation for the functions Tjjj is obtained by subtracting 
Equations (3-19) from Equation (3-14), yielding 
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If the superscript S denotes the short-periodic part of a function such that 




then the preceding equations can be expressed as 
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m- 1,2,...,5) (2-2Sa) 
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Inspection of Equations (3-28) indicates that the functions (i - 1,2,..., 5) 

depend only on quantities of lower order. In the case of the sixtii variable Z , 
the function tjg . also depends on the jth-order function intro<iueed through 

the tei-m Nj . Hence, the function must be determined prior to the function 

n6,j- 

These functions are determined to witiiin aii arbitrary function of the slow var- 
iables, by developing tlie riglit-hand side of Equations (3-28) into multiple 
Fourier series and integrating term by term. More explicitly. 
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The functions . therefore have the form 




(3-30) 


where . is a 2;r periodic function of I with zero mean, i.e., 




(3-31) 


and ia an arbitrary function of Integration dcpentling only on the slowly 
varying mean elements. 

It then follows from averaging Equation (3-30) that 





(3-32) 


This equation is a generalization of the constraints expressed in Equation (3-20) 
and Equation (3-23). Therefore, in order specify the functions . most gen- 
erally, a set of ai^bitrary functions of the slow variables is required. Because 

the function Cj j(a) is an arbitrary function of integration, it can be taken to be 
identically zero, i.e.. 


5 0 


(3-33) 


«tar«by reproducing the coustraiul used to obtain tho form of the avorag.nl oqua- 
tlons of motion givon in liquation (.1-261. Consequently, the validit, of the appli- 
oalloB of ih. constraint eapressed in .ithor Equation (.1-20) or liquation 11-211 
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has been demonstrated. The use of the constraint fiiven \)v Equation 
requires that the j functions be purely short periodic ami or mixed short 
periodic, i.e., 


In summary, a set of functions Hi,] containing only short-{H'rioilic terms can 
be obtained, and the near-identity transformation given by Kijuations is 

OOmpletely specified by the oxpressions 
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3.2.‘l t’omputat ional Procedure 

The determination of the jth-ordcr contribution to the mean cKmiu'iu rates (Kqua- 
tions (3-26) and the functions are mtcrdeivudcnt and must proicinl serially 

on :m order-by-order basis. To illustrate this procedure, the second-order 
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•guattoM ara presaatad more ejqplicittjr. ExpressiOf Eqpatloa <3-26) to aecond 
order yields 


° + 0(«*) (1.1.2 5K,1-35«) 

^ • n » 4 ♦ Mt\ * 4* ^4», 1^(1,!) ♦ # 0(t'> (3-35b) 


Using EquatiOBs (3-13) and the constraints given in Equations (3-22) and (3-23), 
the averaf;ittg operation yields the simplifications 


(n^cvo). • <- i ^ • 0 
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to»pcction of thi* equation indicates that the first-or ter contributions lo tha moan 
element rates. ^ , are independent of the functions rjj ^ However, the 
•eeond-ordcr eontrlbuUons to the mean element rates. require knoudeclKe 

of the fttoctioos TJi, j • Meoce, the cofnputation must proceed as follows; 
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This procedure is followed in extending to higher order the averaged equations 
of motion. 
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3.3 AVERAGED EQUATIONS OF MOTION FOR MULTIPLE PERTURBING 
FUNCTIONS 

The preceding analysis can be extended in a straightforward manner to the case 
of multiple perturbations contributing to each element rate. Examples of such 
cases are: inclusion of more than one spherical harmonic from the nonspherical 
gravitational potential field, multiple third-body perturbations, and 'combinations 
of these effects with atmospheric drag and/or solar radiation pressure. To first 
order in the small parameters, this formulation is identical to summing the first- 
order averaged equations of motion (Equations (3-26)) for each perturbation. 
However, at higher orders, mixed (coupled) terms appear in the averaged eq'.a- 
tions of motion. To illustrate this phenomenon, the case of two perturbing func- 
tions is considered. The corresponding set of differential equations is given by 


« eF; (o.,t) t>G-,(U,0 (i = 1.2 5) (3-41a) 
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The near-identity transformation (Equation (3-3)) is generalized to 
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where the functions ~ j periodic in the mean fast 

variable i . 

The transform of the original system (Equations (3-41)) is assumed to be of the 
form 
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where the functions Bj . , = B. . , (a) depend only on the slowly varying ele- 
ments. Equations (3-43) are a generalization of Equations (3-4) given previously. 

The constraint 1 < j + k is imposed on the lower limits of the double summation 
in Equations (3-42) by the assumption that the difference between the osculating 
and mean elements is, at most, of first order in one the the small parameters, 


I 0.; “ a; I « max [oU'), OCv)] 


(3-44) 


Similarly, the same constraint is imposed on the lower limits of the double sum- 
mation in Equations (3-43) by the assumption that the magnitude of the mean eie- 
ment rates is, at most, of first order in one of the small parameters. 

In Equations (3-42) and (3-43), the upper limit on the summation over i , N , is 
chosen such that all contributions through order O(e^) are retained. Terms with 
increasing powers of u obviously require decreasing powers of e in order to meet the 
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criterioa that only terms of order less than or equal to O(e^) be retained, i.o., 
0(e^ ii O(e^) . Specifically, for a given value of j , the maximum value of 
k , M(j) , is given by the integer part of the expression 


M(p » (N-n ^ 


(3-45) 


and the range of M(j) is Os M(j) s I « for N s j ^ 0 . 

Differentiating Equations (3-42) and substituting Equations (3-43) for the mean 
element rates into the result yields 
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which is a generalization of Equations (3-6). Rearranging Equations (3-46) yields 
the following generalization of Equations (3-7): 
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advanlageous to decompose the above double 


summation as follows: 
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Expanding the perturbing fvmctions and G. in Equation (3-41) as a Taylor 
series and then arranging as a power series in the small paranietor yields 
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etc. The functions g. . ^ are identical to the above, \vith the exception that the 

If J* 

function Fj(a, A) and its partial derivatives are replaced by the function Gj(a, 1) 
and its partial derivatives evaluated at a = a , L = t • 

Expanding the mean motion, n, as a power series in the small parameter yields 
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and so forth. 

Equations (3-49), and (3-51) are substituted into Equations (3-41) and 

terms with like powers of the small parameters are sot equal, tluis obtaining 
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Equations (3-53) and (3-54) are essentially identical in form to the correspond- 
ing equations for a single perturbing function given by Equations (3-14). Conse- 

qiiently, to obtain the complete higher order contributions for the case of two 
perturbing forces, Equations (3-55) representing the coupled terms must be 
added to the equations for each perturbation given by Equations (3-14). Equa- 
tions (3-53), (3-54), and (3-55) are then averaged (essentially as before), yielding 
the averaged equations ot motion; the remainder of the solution then proceeds as 
before. The final results are as follows: 
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(1 < j ^ N-1 ; 1 ^ k < M(j)) 


where the superscript S again denotes the short-period part of the function. Thu 
assumptions used to obtain Equations (3-56) and (3-57) require that 




as in the case for the single perturbing function. 
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The explicit aicraKcd cqualionH of motion to d'^cond order in both small param- 
eters reduce to 




4 

r * 






to • 


»/^ i^N’’ 

\V.‘ as^ 4 


58 a) 


(i ~ 


I 


r 

1. 


^ • n V € 

dt 




E i / , dPfc IS n , a \ 


)*i 


(;i-i)«b) 


^ -^v4> ib \ 

^ ^ + .i£ ^ \ 


and the ordci of cninputation follows as in the single jHTfurbiai; fnru tit»<i eausv. 
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Extension of the discussion to an arbitrary number of perturbing functions is a 
straightforward but tedious exercise. The differential equations take the form 
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(i - 1,2,...,5) (3--9n) 
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where (k - 1.2,. . ,K) are K distinct small parameters (i.e. , ^ > 

^2 = etc.) and H. is the kth perturbing function acting on the ith element 
(i.e., Hj^ j = Fj( a, X) , 2 “ ^i(^» ^ ) » ) • The near-identity transforma- 

tion and the transform of the above differential equations become 
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The perturbing functions are expanded in the form 


etc. Further pursuit of this procedure adds little additional insight except for 
an appreciation of the cumbersome expressions obtained for the final results. 

A less involved approach is presented next, based on the fact that the practical 
application of such theories is almost always limited to at most second order in 
the small parameters (see Section 3.5). As previously shown in the case of two 
perturbing functions, the second-order averaged equations of motion for a single 
pe: 'nrbation are summed and the coupled term is then added. For K perturbing 
functions, the same procedure holds to second order, i.e., K equations of the 
same form as in the single perturbing function case are summed. The coupled 
terms are then evaluated to complete the second-order contributions. The number 
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of coupled terms is simply the distinct number of pairs obtained from K objects 


taken two at i. time, i. e. , 


k(k-i) 



This procedure provides for all contributions from the K perturbing functions to 
the averaged equations of motion through second order in all K small param- 
eters, Yjj . 
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3.4 MODIFICATION OF THE AVERAGING OPERATION FOR RESONANT 
PHENOMENA 

A commensurability of two mean mocions appearing in the dynamical system, e.g. , 
the satellite and third-body mean motions or the satellite mean motion and the 
central body rotation rate, can contribute significantly to the long-period motion 
of the satellite. The generalized method of averaging presented in Sections 3. 2 
and 3.3 is directly applicable to cases involving such resonance phenomena. 

The basic objective in applying the method of averaging to the orbital equations 
of motion is the removal of short-period terms. The averaging procedure de- 
fined by Equation (3-15) removes the high-frequency components of the motion 
for the majority of problems but is not suitable for the treatment of all resonance 
phenomena. In those cases for which resonance phenomena are significant, the 
averaging operation given in Equation (3-15) may have to be modified. The neces- 
sity of this modification depends on the criteria used for selecting short-period 
terms and the characteristics of the perturbing functions. 


3.4.1 Frequency Characteristics Specific to Resonant Phenomena 


The existence of a resonance condition, i.e., a commensurability in the mean 
motions of the fast variables of the perturbed and perturbing bodies, dictates that 
these fast variables cannot be considered mutually independent. An arbitrary 
term in the Fourier sei'ies expansion for the perturbing function takes the general 
furm 


cos()t - Vi +Bj^l,sin(iJ.-kJ.V6j,') (3-63) 


where JL and X' are the fast variables of the perturt 2 d and perturbing botlies and 
and linear combinations of slowly varying angles. 
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The fast variables X and X* are assumed to have the mean motions n and n' , 
respectively. If the ratio of the mean motions is approximately equal to the 
ratio of two integers, i.e. , 


n* 
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N* 


then 


N'n - Hn' 0 

The fast variables thus obey the relationship 

K'l - Nl' • /4. 


(3-64) 
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where the function m ” M(t) is a slowly varying angle which produces only long- 
period effects. 

One of the fast variables can be eliminated from the perturbing function using 
Equation (3-66), resulting in a formulation dependent on only one fast variable 
and an additional slow variable ix{t) . Eliminating the fast variable X' from 
terms of the form given in Equation (3-63) yields arguments of the form 
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Elimination of the fast variable JL in favor of JL' yields trigonometric arguments 
of the same lotm. More specifically, the quantitiCvS N and N' are interchanged, 
X is replaced by X' , and d! is defined by the sum rather than the difference. 

In general, arguments of the form given in Equation (3-67) produce fractional as 
well as integral multiples of the fast variable JL . This is specifically the case 
when kN' is not a multiple of N . An arbitrary decision to consider only integral 
multiples of the fast variable as short period is not practical in this case, partic- 
ularly in view of the desire to maximize the integration step size. For example, 
the case of a close-Earth satellite in a 12:1 i*esonance with the Earth's rotation 
is considered. From Equations (3-64), N = 12 and N' - 1 , and the argument in 
Equation (3-67) can be expressed as 
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This argument will contribute terms containing the fractional arguments 

1/12X, 1/6X, 1/4X, 1/3X, 5/12X, l/2i, 7/12X, 2 3X, 3 4 i, 

and 11/12 £ 

for those values of k which are not multiples of 12. The averaging operation 
defined by Equation (3-15) will not remove terms with these arguments. IXjfining 
terms containing the arguments 2l and I as short period and terms containing 
1/2 i , 11/12X, etc., as long period would restrict the integration step size to 
approximately one-eighth of the satellite revolution {xjriod. To maximize tlie inte- 
gration step size (hopefully to the order of several orbital jxsriods), while retain- 
ing the basic long-period behavior of the dynamical system, all dependence on the 
fast variable should bo eliminated. This requirement is identical in philosophy to 
that imposed in the selection of the averaging operation for nonresonant plienomena 
(Equation (3-15)). 
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3.4.2 The Averaging Operafion for Resonant Phenomena 

When resonant phenomena are included in the equations of motion, the selection 
of an optimal averaging operation is dependent on the form of the perturbing 
function. The resonant contribution is embedded in this function and is isolated 
by the application of the averaging operation to the function. For this discussion, 
the resonant perturbing functions are separated into two categories: embedded 
resonant terms and quasi-isolated resonant terms. These categories are dis- 
tinguished according to whether or not the perturbing functions contribute terms 
with fractional multiples of the fast variable. 

An embedded resonant term contributes fractional multiples of the fast variable. 
Such formulations of the perturbing function are frequently encountered in numer- 
ical averaging applications where the perturbing function is formulated in terms 
of the complete perturbing acceleration (Equation (2-15)). 

The second category of perturbing functions (i.e., quasi-isolated resonant terms) 
contributes only integral multiples of the fast variable. The resonant contribution 
has been partly isolated from tlie complete pertui'bing function such that only inte- 
gral multiples of the fast variable appear. Moi'e specifically, the perturbing func- 
tion is I'estricted such that the integer k in Equation (3-67) takes on only values 
which are multiples of N , i.e. , 


k « pN 


(P = 1,2, 


■ • • 


) 


It is important to note that no restriction has been placed on the integer j in 
Equation (3-67). Since only particular values of j produce tlie I'esonant contri- 
bution, the quasi-isolated resonant term contributes both short-jxM'iod integral 
multiples of the fast variable only) and resonant contributions to the motion. If 
the values of j are restricted appropriately, the resonant term is comp’etely 
isolated from the perturbation function. 
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As an example of a quasi-isolated resonant term, the 12:1 resonance example 
cited previously is again considered. If k is restricted to multiples of N , i.e., 
multiples of 12, then all fractional multiples of the fast variable are eliminated. 
These terms correspond to the tesseral harmonic terms in the geopotential of 
order 12. In this case, any geopotential term of order 12 would be a quasi- 
isolated resonant term. The specific resonant term, which will be isolated by 
the application of the averaging operation, corresponds in this case to the value 
of j where j = 1 . 

3. 4. 2.1 The averaging Operation for Embedded Resonant Terms 

In the case of embedded resonant terms, fractional multiples of the fast variable 
appear in the perturbing function. In view of the form of the argument given in 
Equation (3-67), all dependence on the fast variable Ji can be removed by defin- 
ing the averaging operation to be the definite integral over the angle cr = A/N on 
the interval 0 < cr < 27T . Expressing a function of two fast variables denoted by 
H in terms of the fast variable a and the slow variable p yields 

* H' (oL,£, yu.) = K*(CL,0',yu.) 


The average of the function H*(a, a, p ) is defined as 





,0-,^) do- 


(3-68) 
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The averaging definition can be expressed explicitly in terms of the fast variable 
jt. IfO<(T^2;r, then 0 < i 1 2 n’N and 
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(3-69) 


Therefore, in the case of an embedded i '.onant term, the definition of the aver- 
aging operation should be specified as 



^3-70) 



This definition has been used by Si'hubart (Reference 43) for performing a numer- 
ical investigation of the Hilda group of minor planets which exhibit a 3:2 conimen- 
surability with Jupiter. Also, ikmson and Williams (Reference 44) used tlie same 
definition in their numerical invest' gat ion of resonances in the Neptime-Pluto 
system. 

It should be noted that tlie above averaging operation removes only those terms 
with periods of 2n’N or less. It does not remove any contributions to the motion 
caused by the resonance, since the fundamental |Xjriod in the motion caused by 
the resonance is contributed by the angular variable p ;md is given by 

9l7t 

"n'm - Nn' 
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3.4. 2.2 The Averaging Operation for Quasi-Isolated Resonant Terms 

Since only integral multiples of the fast variable, X , appear in the case of quasi- 
isolated resonant terms, the averaging operation given in Equation (3-15) is 
applicable. It is repeated here for convenience: 



(3-71) 


where H* denotes a quasi-isolated resonant term. 

The distinction in the averaging operations given in Equations (3-70) and (3-71) 
has an important implication for numerical averaging theories whore the aver- 
aging is performed using a numerical quadrature. The perturbation model must 
be evaluated at each abscissa in the quadrature interval (usually between 12 and 
.)6 points per interval). Numerieally averaging an embedded resonant term 
requii’es N times as many force evaluations as the numerical averaging of a 
quasi-isolated resonant term for a total of Ixitwcen 12N and 3GN force evalua- 
tions. Therefore, in the application of the numerical averaging methods, the 
perturbation models should be restricted to the quasi-isolated resonant terms 
whenever possible. 

The spherical harmonic expansion representing the nonspherical gravitational 
potential is well suited for obtaining by inspection the ciuasi-isolatcd resonant 
terms. The commensurability is directly related to the order of (liose terms 
which contribute to the resonance. Such is not the case for the closed-form, 
third-body perturbing acu’oieration or even foi’ the standard expansion in 



Legendre polynomials for the third-body disturbing function. The resonance 
contributions remain embedded in these particular forms. However, the third- 
body disturbing function can be expanded in spherical harmonics using the asso- 
ciated Legendre polynomials (Reference 18). The quasi-isolated resonant terms 
are then immediately obvious as in the case of the nonspherical gravitational 
potential. 
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3.5 THE APPLICATION OF HIGHER ORDER AVERAGING THEORIES 

The implementation of a jth-order theory requires the explicit determination of 
the near-identity transformation through order j-1. Consequently, higher order 
averaging theories are significantly more complex than the first-order theory. 
Examination of the second-order averaged equations of motion indicates that a 
first-order theory should suffice for all cases where the amplitude of the first- 
order short-period variations in the osculating elements are small, either abso- 
lutely or relative to the amplitude of the long-period variations in the mean 
elements. 

In cases where a second-order theory is needed, it should be applied selectively 
to those terms producing the largest short-period perturbations, e.g. , the 
oblateness (Jg) term in the zonal harmonic expansion or the first few terms in 
the expansion of the third-body disturbing function. Such restrictions are usually 
justified on physical grounds and by the practical considerations of implementing 
a higher order theory. For those cases where such restrictions cannot be 
justified on physical grounds, an alternate formulation of a problem, e. g. , a 
restricted three-body problein, should be considered. 

3.5.1 The Significance of Second-Order Terms 

Two questions are of particular interest concerning the possible significance of 
second-order terms in the averaged equations of motion: 

• How do the solutions of the first-order equations and second-order 
equations differ with time ? 

• What are sufficient conditions such that second-order terms can be 
neglected over the time interval 0 ^ t < T ? 

A precise answer to the first question is impossible without generating the actual 
solutions: however, a qualitative estimate of this behavior is possible. The 
answer to the second question is provided by inspection of the second-order aver- 
agod equations of motion. 
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3. 5. 1.1 A Qu.ilitative Comparison of the First- and Second-Order Theories^ 

The quantity [ a'(t), J'(t)] is defined to be the solution of the following system 
of second-order averaged equations: 





0 = 1,2,...,5) (3-72a) 


« n(tXi) e.At(5') + 


t3-72b> 


Similarly, [a*(t), X*(t)] designates the solution of the system of first-order 
averaged equations 





(i - 1,2,...,;")) (3-73a) 


iL ‘ 

di 


nCcL*J 


(3-73b) 


The difference of ihe solutions is designate;! as 


r,CO » O-UO - 


(i 1,2,...,.")) 


vj.i) -- I ’ Ct ) - i *(0 


t;5-7-aii 


This discussion follows closely that given by \\ . T. Kvner in a s-aries of lectures 
on the topic of nonline;M' resonance tsee Ueference Si. 
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A set of differential equations representing these differences is given by 


£ [A;,i(a') - A; Jto] 


(i -• 1,2,..., 5) (3-75a) 


n(3.'0 - n(aj) fe [Afc,i(t') - A^.^d*)] (3-75b) 


It follows that 




(3-76a) 


and, similarly, 


jr^jCOI < f I ndi) - I dt' 

-t 

+ I (a.') - d*) I dt' •*• £*“ f I Afc^iCa.*^ I dt' 

^ A » A 


,3-76b) 


The functions j and n are assumed to satisfy the Lipschitz condition 


I A\,id') - A, < Lili'-t*! 5 L-,|?(t'l| (i = l,2,...,fi) (3-77a) 


ORIGINAL PAGE IS 
OF POOR QUALITY 




f r 

♦ • 


-♦ ^ , 




’ I 
: I 

• S: 


I 

I 

r 

• ( 

1 : 

1 J 
1 

'1 


I 


I - r\('oJ\) I i L' |L'i- (ttl L' |a'- 0.*| = L | r(t) | c>,-77b) 


on I interval 0 < t < T , where U and L' are ;x)sitive condtants and wlicre the 
vector ‘r consists of the components I’j (where i = 1,2,..., 5). It is sufficient 
that the partial derivatives of the functions /V; , and n exist and arc l»imdcd on 
the interval 0 < t < T for Liquations (.‘1-77) to be satisfied. It is also assumed 
that the absolute value of the second-order function A <, is bur . • iI>ovc 

on the interval U < t < T , i.e. , 

|A;^i| < Mi for 0< t < T 


Substituting tiquatiois (;i-77) into Liquations (;l-7G) yields the inequalilies 




\T;ii)\ < cL; j |fCn| at' C^Mit 

C a 

|r^(x')| < €L^,/ |fa')|di' 4 i‘l |fa')|dt‘ 4£M(,t 


(1 1,2,...,.) );i-7sa) 


(:i- 7sl. 
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Then, summing EquK.ions (3-78a) over i yields the inequality 

^ t 

|-?a)| < y |r,Ct)| < €L/ IrttOiat' 4- eU 

1*1 


n also follows that 


I r^a') I i c iiv £') L I (f 1 aV 


£^Lt 


Using the generalized Gronwall inequality,^ it is easily shown that 


x: ^ 

jfa)l < J' fe^Ldr 

0 

= f L eApTfeUi-x')] dr 

'^0 

= £ [cKp(6.Lt') - l] ^ £^ Lt e)kp (feLt) 


^The Generalized Gronwall Inequality (Reference 45) 


If the following four conditions are met: 

(1) > ^t), </>(t), and u(t) are defined on the interval Iq < t < T 

(2) A.(t) is greater than or equal to zero and is summable 

(3) <^(t) and u(t) are absolutely continuous 

(4) the following inequality is satisfied 


u,(0 ^ J Mz) a(x) dr + 


(tQ< l<tj) 


1 

U.CO 6 Mr) dz^ J do) 

dz 
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Substituting the miaimum of the upper bounds for |'r(t) | , i.e., 


i. 6 [exp(€Lt) - i j 


into the inequality for | r (t) | yields 


i (lv€')6Lj[cxp(€Lt')-l]dt» 


•I- 6^U 


(3-81) 


= (l>t) [ftXp(eLO - 1 - €Lt] + e^U (3-82) 


s €Lt (€ -*• feLt) exp (eUt) 


It is noted that this last result is not in agreement with that obtained by Kyner, i.e. , 


< eUt [3LeAp(6Lt) + 6-t] 


(3-83) 


In summary, the difference of the first- and second-order solutions is bounded by 
the functions 


I U* - a.* I i e* Lt e.Ap (glO 


(3-84a) 


exp(€lO 


; .1 

1:1 • 
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If the time t is restricted such that eLt«l (i.e., Os t s T« (eL) ^), then 
generally the order of magnitude estimate of the divergence between the first- 
and second-order theories is given by 


-o(gH) (for 0^ t i T«(€L)“^) 


(3-85a) 
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(3-85b) 


|i'- i*| -o(dH*) 


(for L“^< t sT«(eL)"^) (3-85c) 


The above error estimates can be mapped back into the osculating elements using 
the near-identity transformation. Only first-order terms are assumed, since 
only first-order terms are required for the second-order averaged equations of 
motion. Evaluating the near-identity transformation 

i » O.^ 0(fe*) 

X » +OCe*) 

with the elements obtained from the first- and second-order solutions and taking 
the absolute value of the difference yields the inequalities 

4 p-sea) 

|i'-x*l i ir-I'l (3-86b) 
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where is a vector with the components y^. ^ (i -• 1, 2, . . , , 5). 

Since the constant L can be chosen to satisfy the Lipschitz conditions 

|ntS',l')-^(S*I«)l 6 LlfW|»L|rtCW| (3-87a) 

* L|f(tt| t LlftCOl (3-87W 

Equations (3-86) can be simplified to give 

S (1+60 |‘r(t)| + feLj r^Ct) j (3-88a) 

I A* -A* I ^ feLlF (t)| + (l+eL)|r^(Ol (3-88b) 

Substituting the upper bounds for | r(t) j and | rg(t) | into Equations (3-88) yields 
the inequalities 

I X' - a* 1 <1 Lt ( 1 + 3L€L ^r 6 L^i') CKp (felt) (3-89a) 

I a' - X* I i 6^ Lt (l+ a£ L + e.l\ + L\) eJtp (tit) (3-89b) 

which yields the following qualitative estimates for the osculating elements: 

lt'-a*| ~ 0(4»t) (for 0 ^ t < T « ( € L)"^) (3-90a) 
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(for 0 s t s L“^) 


(for 0 s t s T«(cL"^) 


(3-90b) 


(3-90c) 


when the restriction €Lt«l is imposed. Thus, the qualitative behavior of the 
oscu'ating elements is, in general, the same as that of the mean elements (Equa- 
tions (3-85)). 

In summary, for arbitrarily small € , the difference in the first- and second- 
order theories is arbitrarily small. For a given e , the difference in the theories 
will be sufficiently small for some interval of time 0 s t sT , where T ~ 0(e"^). 

3. 5. 1. 2 Sufficient Conditions for Neglecting Second-Order Terms 
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The second-order averaged equations of motion are given by 

^ A-, gid) + OCe^) (i=l,2 5) (3-91a) 

• n(aO + £A6^^(a) ^ O(fe^) (3-9ib) 


where 



(i 1,2,..., 6) (3— 92a) 


(i = l,2 5) (3-92b) 
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(3-92c) 


Inspection of tJu? second-order averaged equations of motion indicates that, for 
the limiting case in which the first-order short-period variations of the osculat- 
ing elements are identically equal to zero, the second-order contributions to the 
mean element rates vanish identically. Similarly, if the amplit- 'cs of the first- 
order short-period variations are small in magnitude, the second-order contri- 
bution to the mean element rates will be small, provided that the short-periodic 
part of the function dFj/Sajj is not large. Finally, inspection of the second- 
order equations indicates that the effect of nonzei'o second-order terms will be 
most significant when the first-order contribution to the mean element rates is 
very small or zero. Consequently, the inadequacy of a first-order theory will 
be most apparent when the element history approaches a local maximum or 
minimum value. 

Before further discussion, the following relation will be demonstrated: 

dFid.i) dA-n(a) 

45-v “ daw. div 


w’here ( * ) indicates d( )/dt. (Since extension of the following disc ussion to 
the case i = 6 is straightforward, it is not presented. ) 

Substituting the relation 


F- CoL,0 * F;(t J) -V OU) 


(i = 1,2,. ..,5) (3-94) 
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into the high-precision equation (from Equation (3-2)) 


original. 
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page 

qUAUT^' 


dt 


6 F-, ( a.,0 


(i = 1,2,...,5) (3-95) 


yields the result 





(i = 1,2,...,5) (3-96) 


Differentiating with respect to time the near-identity transformation from Equa- 
tion (3-3) 



^ ilia 

dt * dt 


(1,1) i. OCe^) 


(3-97) 


and substituting into the result the expansion of the mean element rate from Equa- 
tion (3-4), i.e. , 


"5T 


6 Ai,iCt> + OCe») 


(3-98) 


yields the relation 

^ ^ + OU») (3-99) 
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Comparison of Equations (3-96) and (3-99) yields tlie result 





at 


(i == 1,2,...,5) (2-100) 


and Equation (3-93) follows immediately. As a result, the second-order contri- 
bution to the mean element rates reduces to 



The requirement that the magnitude of the short-period part of the function 
3Fj, 5aj^ is not large then reduces to the requirement that tlie magnitude of tlie 
function 1 is not large. It seems reasonable to ex{X'ct that, if tlie 
function j has a small absolute variation and tliore arc few local extrema 
over the interval i-orresponiling to one satellite period, then the first time der- 
ivative of the function should not be large. This assumption shouUi also extend 
to the partial derivatives. 
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A somewhat more formal criterion for neglecting second-order terms requires 
simply that the integrated effect of the second-order term over the interval 
0 < 1 1 T be less than some specified tolerance 6 , i.e.. 




(i = 1,2,...,6) (3-103) 


fid 


or, more specifically (in view of Equation (3-101)), 




(3-104) 


Clearly, the integral of the second-order contribution can be bounded as follows: 


|a,,3lC^) dt 


and it follows from Equation (3-101) that 


(3-105) 


i ‘ I: I <^... 
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(i = l,2 5) (3-106) 
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where and iTj designate the maximum variations of Uie functions j ami 
, respectively, i.e.. 
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For the case of the fast variable (i.e. , where i - G), 


(3-107a) 


(i ^ (3-107b) 


^6,^ ■ 
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(3-108) 


and, therefore. 




( 3 - 109 ) 
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It follows that 


t t 

6^A,,(l^dl < . tM,t (i = 1,2, .. ,,G) (3-110) 
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Thus, the second-order term may be neglected when e^M. t < 6 (i - 1,2,..., 6), 

o 

that is, over the interval 0 < t < T , where T = 6/(e Mj) . Therefore, the time 
interval over which a first-order theory is valid depends inversely on ihe magni- 
tude of the first-order short-period variations in the osculating elem«'. 

A relative criterion for neglecting the second -order terms provides .'-r n ..ttle 
more insight in practical applications. Essentially, it is required that the inte- 
grated second-order contribution be negligible when compared with the integrated 
first-order contribution over the interval 0 < t < T . Specifically, the condition 
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(0 t < T1 


(3-111) 


is to be satisfied.^ As before, the integral of the second-order term is easily 


bounded by 


J A',,1 (H S M, t 


Also, if die following definition is made 


(3-112) 
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(i = 1,2,...,6) (3-113) 


The coupling between the second-order and first-order contributions is assumed 
to be negligible. This argument is valid only for the bounded periotlic elements 
or very slowly gi'owing secular elements, since the rapid first-order secular 
growth of the fast variable would satisfy the inequality even for lai'ge second- 
order contributions. This critcri.>r. is really useful as a negative criterion spec- 
ifying when second-order terms an lefinitely neces ary rather than when they 
can be neglected. 
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If (Equation (3-107b)) is replaced by the order of magnitude estimate 



^n..i 



(3-116) 


where A denotes the maximum variation of the element over the interval 
0 < t' < t and p! is defined to be an upper bound of the time derivative of the 
short-period variation j , i.e., 

lauU/o; 

Then it follows that 
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and the second-order terms can be neplected when 
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Thus, the second-order terms can be neglected over the interval 0 <' t < T , where 


T « 


Z— I 


(3-120) 
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The smaller the ratio of the short-period variation to the first-order long-period 
variation, the greater the interval over which a first-order theory is valid. Hov/ 
small these ratios must be depends on the time interval over which the first- 
order theory is to be valid. The answer to this question can be provided only by 
a thorough investigation for each dynamical system. However, an upper bound 
of at most a few percent would be a likely guess for retaining a period of validity 
of a few years. 

On the other hand, a first-order theory is clearly inadequate when the amplitude 
of the short-period variations is 20 to 30 percent of the long-period variations. 
The author has investigated the case of a near-circular satellite (IMP-J) in 
2:1 resonance with the Moon. The amplitude of the short-period variations was 
approximately 30 percent of the magnitude of the long-period variation caused 
by the resonance. The first-order averaging theory produced poor results in 
the acighborhood of a local extremum of the semimajor axis history, an indi- 
cation of significant second-order contributions to the motion. 

3.5.2 Application of a Restricted Second-Order Theory of Averaging 

The application of a second-order averaging theory to all perturbations would 
compromise the advantage of the low computational cost, which is characteristic 
of the first-order theory. However, the application of a second-order averaging 
theory, restricted to selected perturbations, may yield more accurate results 
where the application of a first-order theory is marginal, or it may extend the 
time interval over which the first-order theory is valid with a minimal increase 
in cost. 
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3.5.2, 1 Nonspherical Gravitational Perturbation 

The spherical harmonic expansion representing the potential of the nonspherical 
gravitational field of tlie central body (Earth, Moon) contains the small param- 
eters 



where a^ designates the equatorial radius of the central body, the quantity a 
designates the semimajor axis of the satellite orbit, and the coefficients 
and Sjj are observed quantities . 

The zonal harmonic coefficients J are defined by 

n 



These small parameters are obviously bounded above by the numerical coeffi- 
cients and . Since Jg » 0(10 '^) and since all other coefficients 

are of the order of for n > 2 , the oblateness term, J 2 , in the geopotential 
might seem to be a logical candidate for the application of a second -order aver- 
aging procedure. In fact, a consistency argument is often made that second- 
order oblateness contributions should be included if any other terms in the 
spherical harmonic expansion for the geopotential model ax'e also included. 
According to the previous discussion, this is not necessarily the case since the 
second-order contributions depend on the first-order shnrt-period variations of 
the osculating elements and their time derivatives and not on the first-order con- 
tribution to the long-period motion. However, it is reasonable to expect that if 
second-order terms are necessary, the J 2 contribution would strongly dominate 
over the other harmonics. Consequently, any second-order theory for the non- 
spherical gravitational field could be limited, in most cases, to the Jg oblateness 
contribution. 
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3.5. 2.2 Third-Body Perturbation 

The case for the third-body perturbation is not generally as simple. The rele- 
vant small parameters are the nth power of the parallax factor, i.e. , 

(n = 2, 3, . . . ) 

where a and a' are the semimajor axes of the satellite orbit and third-body 
orbit, respectively. (It is tacitly assumed that the disturbing third body is 
an exterior perturbation, i.e., a < a' . If, however, a > a' , the expaJision 
proceeds in powers of the inverse of the above parameter. ) 

The upper bound of this set of small parameters is unity in contrast to the upper 
bound for the small parameters in the nonspherical gravitational model which is 
of the order 0(10“'^). Clearly, for high-altitude satjllites, the small parameters 
are not really very small except for the large values of n. Physically, as the 
parallax factor grows towrad unity, the third body produces stronger disturbances 
(botli short- and long-period) in tlie satellite motion. These larger disturbances 
require a more complex model which is manifested by a greater number of terms 
in the disturbing function expansion. 

The recursive formulation of tlie disturbing function presented in Volume II of 
this report ciin, in principle, be used to produce expansions to any arbitrary 
order. However, high-order expansions can nroduc'e increased computational 
cost, unavoidable numerical round-off and truncation errors, and, possibly, 
errors due to unstable recursion formulas. Also, the first-order averaged 
equations of motion can be formulated in terms of the perturbing acceleration 
using tlie Gaussian formulation (Kquation (2-15)) to avoid entirely the problem 
of slow convergence of the disturbing function expansion. 

The slow convergence of the disturbing function is of far greater significimee 
to t:e application of the method of averaging itself. The strong short-period 
disturbances can no longer be neglected in formulating the averaged equations 
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of motion; a higher-order theory becomes necessary. The stronger the short- 
period disturbances, the greater the number of terms in the disturbing function 
expansion which must be developed to second or higher order in the application 
of the method of averaging. 

An additional complexity is presented because the small parameters are not 
entirely independent, specifically, 

m 

^ mm 

Under these constraints, an order of magnitude argument would indicate that if 
the term containing the small parameter eg produces significant short-period 

9 o 

contributions to the motion, then the terms containing the powers 

and the products Cg * ^2 * ^4 contribute significant short-period 

variations and, therefore, f'ould be required in the averaged equations of motion. 

On the other hand, the numerical coefficients in the higher order terms of tlie 
averaged equations of motion may render their contributions less significant 
than tlie order of magnitude argument would indicate. If this were generally 
true, the slow convergence of the disturbing function would have a less severe 
impact on the order required in the application of the method of averaging, and 
development of the first few terms of the disturbing function expansion to second 
order might considerably extend the range of the parallax factor win l e the aver- 
aged equations of motion are valid. At the very least, it would extend the appli- 
cation of the method of averaging to those cases where a first-order application 
was only marginal and/or extend the interval of validity to tens of years^ in those 
cases where a first-order theory already provides adequate results over a few 
years. 


Such very long predictions may be necessary to meet futux'e mission analysis re- 
quirements, e.g. , permanent space station missions, solar power satellites, etc. 
High-precision techniques are not well suited for such investigations. 
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Furthermore, even though rough qualitative estimates of the behavior of the aver- 
aged equations is known (see Section 3.5.1), very little is known quantitatively of 
the time intervals over which a first-order application of the method of averaging 
is valid. The basic difficulty in ascertaining an accurate estimate of this time 
interval is obtaining a suitable reference solution with which to gauge the aver- 
aging theory. The standard approach has been to compare the results obtained 
from the averaging theory with a high-precision reference solution directly or 
with mean elements obtained in some manner from the high-precision reference. 
Beyond the time intervals over which high-precision technique^ a.e valid, the 
only reference is that obtained by direct observation.^ 

It is assumed that these mean elements continue to provide an accurate picture 
of the long-period and secular motion of the dynamical system for several years 
or more. Of course, the exact length of this time interval depends on the mag- 
nitude of the short-period variations in the osculating elements as shown in Sec- 
tion 3. 5. 1. 2 of this document. Eventually, the prediction from the first-order 
theory will gradually diverge from the real solution. Without some comparison, 
this divergence will probably become apparent only after it has reached extreme 
proportions. A comparison with a complete first-order averaging theory aug- 
mented to include the dominant second-order contributions from the oblateness 
and third-body perturbations would provide some insight into the period of validity 
of the first-order theory in addition to extending it. This approach. In essence, 
computes the value of the integral in Equation 3-103). 

In summai'.v, a first-order application of the method of averaging is adequate for 
several years in those cases where the short-period variations in the osculating 
elements are either absolutely small or small relative to the long-period varia- 
tions of the first-order mean elements. The implementation of a complete second 


This situation is not peculiar to the averaged orbit generator but also affects 
the high-precision generator, since the only reference is provided by observa- 
tions. 


3-73 








or higher order averaging theory would reduce the computational advantages 
characteristic of the first-order implementation. However, a second-order 
implementation, restricted to only the dominant perturbations, would extend 
the application of the method of averaging to a wider class of problems while 
minimizing the additional computational cost, and it would also provide some 
estimate of the time interval over which first-order averaging is adequate. 
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SECTION 4 - FIRST-ORDER SHORT-PKRlOn CONTRIBUTIONS 
TO THE OSCULATING ELEMENTS 


For many plications, the solution to Equations <3-2) (the true instantaneous 
or osculating elements) is desired. Several techniques have been developed for 
the solution of these equations (e.g. , Cowell, Encke, etc. — see Reference 29); 
however, these high-precision techniques share the characteristic of high com- 
putational cost. To reduce this cost, the averaged equations of motion wore 
developed, which provide moan elements for the dynamical system. 

In addition to the mean trajectory, the method of averaging provides (in principle) 
a way to compute a jth-order approximation to the osculating elements from the 
mean elements. First-order, and possibly second-order, approximations to the 
osculating elements are sufficiently accurate for most applications. The com- 
putational complexity of these approximations increases tremendously with the 
order of the small parameter. 

The effectiveness of representing osculating elements by applying a first-order 
short-period variation to mean elements has been demonstrated by Lutsky and 
Uphoff (Reference 5). It might appear that such a proceduix; would vitiate the 
computational advantages associated with the methotl of averaging, and it has 
already been demonstrated that mean elements are sufficiently accurate for many 
applications (References 4 and 9). However, for some applications, e.g., defin- 
itive orbit determination procedures, the additional accuracy provided by the 
first-order short-period variations might be necessary. 

Based on the following discussion, it appears Utat tlm cost of evaluating the first- 
order short-perloil variations using an analytical formulation^ would be no more 
costly than a single averaged derivative evaluation. This estimate is based on tlie 
assumption that the evaluation of first-order short-period variations is performed 
independently of the derivative evaluation. 

^The cost of evaluating the first-order short-jMJriod variations by a numerical 
technique can be estimated by reviewing the method presentetl in Reference 5. 
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As will be shown in Volume II, the mathematical formalism for the mean element 
rates is also common to the first-order short-period variations. Consequently, 
it is estimated that, if proper advantage is taken of thlG commonality, the cost of 
evaluating the analytical formulation of the first-order short-period variations 
could be reduced to possibly 20 percent (or even less) of the cost of a derivative 
evaluation. For those environments where the utmost computational efficiency 
is required, these variations should be applied only at judiciously selected points 
along and/or at the end of the trajectory for applications with hi^-accuracy 
requirements. 

An equally important application for such approximations to the osculating elements 
Is the conversion of osculating elements to mean elements. An osculating-to-mean 
element conversion can be developed by inverting the equations which specify the 
mean-to-osculating element transformation. 

The mean elements describing the long-period variations in the trajectory are 
only as accurate as the initial mean elements and, hence, only as accurate as 
the osculating-to-mean conversion. Existing conversion procedures are strictly 
numerical (except for the Brouwer theory, which is limited to the low-order zonal 
perturbations) and are based on quadratures or costly differential correction 
procedures which require a high-precision orbit generator. Therefore, either 
the initial conditions must be predetermined or the softwai'e system must have 
access to a high-precision orbit generator as well as to the averaged orbit gen- 
erator.^ In addition, implementation of the short-period corrections appears to 
require no additional theory beyond that necessary for the averaged equations of 
motion. 

This section presents a discussion of the fii'st-order short-period variations of 
the osculating elements and their application to both osculating-to-mean and mean- 
to-osculating element conversions. This discussion is developed in the context of 


If nonconservative perturbing forces, e.g. , drag, etc., are present, there is 
no recourse (at present) to the numerical osculating-to-mean conversions. 
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the Ln^jninfte Planetary t}quations (Kquatlons {2-:U)). A dissuasion of the first- 
order short-period variations in the rontext of tlw (laussi^in \’ariaiion o» Param- 
eters (\’OP) equations Jind the numeriral averaging approach I'an U' found in 
Uefoi'oiu'o 5. 
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4.1 MEAN-TO-OSCULATING ELEMENT CONVERSION 

The near-identity transformation (Equation (3-3)) establishes the relation between 
mean elements and osculating elements. A general expression for the jth-order 
term in this transformation is given in Equation (3-29). Evaluation of this expres- 
sion for hi^er orders is quite complicated if not prohibitive. However, evaluation 
of the first-order term is manageable. (This term also appears in the formulation 
of the second-order averaged equations of motion (Equations (3-39).) 

Expressing the near-identity transformation to first order in the small parameter 
yields 


0i\ =■ Olj + 1»2,...,5) (4-la) 


where 




(4- lb) 


feTJ. ^ J'e (i = l,2....,5) (4-2a) 


(4-2b) 


and Fj denotes the short-periodic part of the perturbing function, i.e.. 


F? * F-. - <Fi)j 
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If the functions 1),^ ^ (where i = 1 through 6) can be evaluated, then Equations (4-1) 
provide a first-order mean-to-osculating element conversion. 

Using the Lagrange Planetary Equations (Equation (2-31)), it follows from Equa- 
tion (4-3) that 




(4-4) 


Because the nonzero Poisson Brackets are independent of the fast variable (see 
Appendix A), 


<■ 


/X .. V mt,h\ _ , V /aR(t.l)' 


(4-5) 


Consequently, Equation (4-4) can be expressed as 


£Fi(t,l).-> (Ui.ZiO 


’L ■ v^/i J 


(4-6) 


Since 
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Equation (4-6) can be simplified to read 


where 


- (r)j 


Substituting Equations (4-8) into Equations (4-2) and simplifying yields 
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(4-lOa) 
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Since the disturbing function R is assumed to be appropriately continuous and 
differentiable, 



(4-11) 


If the short-periodic function S( a, i ) is defined as 


set,, I) » 


dl 


(4-12) 

> 


then Equations (4-10) take the form 






dS 

6 ^ 


(i -- 1,2,...,5) (4-13a) 
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(4- 13b) 


Equations (4-13) are almost identical in form to the general form of the Lagrange 
Planetary Equations (Equations (2-28)), with the exception of the reciprocal aver- 
age mean motion factor and the second term in the equation for the short-|xjriod 




^ Fxnressin Equations (4-13) explicitly 

variation of the fast variable, i * 

L 1 ^ \ \ rp<»ul'a in the following: 

in equinoctial elements ( a, h, k, p, q, X ) 


azc ds 


* nk ~S^ 


(4- 14a) 


S B /dS 

^ I ^ 


arlABX^^P 


(4- 14b) 
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B ^ 

?\A \dh 


hC /- ^ . Q ^ 
anAB V ^ 


(4-l4c) 


pC 1 

L dS 

W -T 


InAB 

V 6V\ 

dk 

IC^ 

55 


4rlA6 




(4-14d) 
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- ^ B /- 6S , ds\ 

?>A 65, nA(UB) 6W ‘6k/ 

(4-14H 

»RAB dp ^ d|y RA S 


where 

S = S(S,A) 

A » n 5.^ 

B - VT-h’^-k^ 

C « 

X * the retrograde factor 


The iadefiaite integral of Equation (4.Ma) yields the last term in Equat.on (4-14f) 

These equalioos ea„ also be expressed In U,r„,s of (he diree.ion .-osines ,5. S.t) 
through Equations (2-d5) and (2-36). 

ExpUelt compufaelon of s for «,„spherU al gi avilallonal pern.rh.ng fu,K ..on 
and fer Ihe thiid-boil.v perturbing funetlon is d.seussed in Volume II of this report 
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4.2 OSCULATING-TO-MEAN ELEMENT CONVERSION 

An osculating-to-mean element conversion is immediately obtained by inverting 
Equaticrs (4-1). These equations are identical in form to Kepler's equation and 
can be numerically inverted, i.e. , solved for the mean elements, by the same 
techniques. These techniques require an iterative scheme, since these Kepler- 
type equations are transcendental. 

Expressions for the mean elements are obtained by writing Equations (4-1) in the 
form 


0.-, = £L-, - 
and 

I • i - 67^4 ; ci,i) (4-15b) 


An a priori estimate of the mean elements will permit evaluation of the right-hand 
sides of Equations (4-15). This, in turn, permits a computed approximation to 
the mean elements. These approximate mean elements are used to reevaluate 
the right-hand sides of the equations and to compute a new approximation to the 
mean elements. The kth approximation to the mean elements is expressed simply 
as 


5-i.k ' 5) (. -Ca) 

^k * ^ ^k-L^ (4-16b) 
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Such a procedure should converge within two or three iterations, provided a good 
a priori estimate is used. 

A good estimate of the initial mean elements is provided by the osculating ele« 
ments. The osculating elements differ from the true mean elements by order 
e and, hence, introduce an error of only second order (1. e. , 0(c^)) when used 
to evaluate the right-hand sides of Equations (4-16). 

It should be noted tha* Equations (4-15) are of the same form as those given by 
Brouwer (Reference 44) and for transforming from the Brouwer primed element 
set (containing the long-period and secular motion) to the Brouwer unprimed 
element set (a first-order approximation to the osculating elements). Equa- 
tions (4-lS) take on the more familiar form of Brouwer's formulas when the 
expression for the functions , given in Equations (4-14), are introduced. 
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APPENDIX A - THE EQUINOCTIAL ELEMENT SET 
AND REFERENCE SYSTEM 

A.l DEFINITION OF THE EQUINOCTIAL ELEMENT SET 

The equinoctial elements defined in terms of the Keplerian or classical elements 
are given by 


Ou a Ou 

V\ a e sin (a) fin) 
k a e cos(tO f m) 

(A-l) 

p a tan^Ci/a.) sinO. 
a "tar?" 0|a') co£»fl 

A » i. + tdfin 


where I is the retrograde factor and assumes the values 

I a 1 for 0 < i < 7 t/ 2 
I a for 7 t/ 2 < i < ;r 

If I = 1, the resulting element set is referred to as the direct equinoctial elements 
and for I = -1 the retrograde equinoctial elements are obtained. The direct equi- 
noctial element set produces a singularity in the Variation of Parameters (VOP) 
equations for the inclination value i = ;r and the retrograde element set produces 
a singularity for the inclination value i = 0 . Hence, both element sets are re- 
quired if the possibility of a singularity in the VOP equations is to be avoided. 
Since the inclination value i = tt is seldom encountered, the direct elements will 
suffice for the vast majority of applications. 

Defining the value of the retrograde factor based on the cut-off value i = tt/ 2 is 
quite arbitrary, and there is no compelling reason to change from direct to 
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retrograde elements (or vice versa) in the middle of a numerical integration 
simply because the value of the inclination passed through this arbitrary cut-off 
value. On the contrary, this cut-off value is intended only as a guideline for 
choosing, at the initiation of the integration procedure, the element set to be 
used. 

In Equations (A-1), the elements h and k are the components in the appropriate 
(direct or retrograde) orbital frame of the eccentric vector, with magnitude e , 
directed toward the periapse. The elements p and q can be considered as 
the components of a vector with magnitude tan(i/2) directed toward the ascend- 
ing node. The element X is the mean longitude. 

Equations (A-1) are easily inverted to provide the transformation from the equi- 
noctial to the classical elements, i.e., 


CL ■ Cl 


£ * VhHl? 


t 

\ * 


(i) » 


arcos 


arctavi 


I arctan ( 

\V 



Cl * ftrctan 

I * X - - in 



(A-2) 
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A. 2 THE EQUINOCTIAL REFERENCE SYSTEMS 


The equinoctial reference frames (direct and retrograde), designated by the 
orthogonal triad (f^ g, w) , are right-hand systems and use the satellite orbit 
plane as the fundamental plane of reference. The unit vector t is directed 
toward a point in the satellite orbit displaced from the ascending node through 
the angle -n for the direct system and through the angle R for the retrograde 
system. The unit vector w points toward the north equinoctial pole and is 
identically the unit angular momentum vector. The vector g is directed toward 

/V 

a point in the orbital plane 90 degrees in advance of the unit vector f and can 
be expressed as 

^ X W X T 

The relationship between the equinoctial I'eference systems and an arbitrai'y 
right-hand reference system, e.g. , the equatorial system, is shown in Figures 
A-1 and A-2. Clearly, in both the direct and retrograde cases, a series of 
throe rotations is required to make the arbitrary reference system coincide 
with each of the equinoctial reference systems. More sjiecifically, a positive 
rotation about the z axis through the angle S2 points the x axis towarti the as- 
cending node. A positive rotation about this new x axis tlirough the inclination 
angle, i, rotates the x,y plane into the f, g plane. Finally, for tlie direct 
case, a rotation about the curi'ent z axis (coincident with tlie w vector) through 
the angle -fl points the x axis along the f vector. For the I'otrograde case, 
this last roatation about the z axis is performed through the angle R to align 
the X a.xis with the f vector of the retrograde system. This series of rotations 
pi'ovides the transformation of the coordinates of any point (e.g., satellite pos- 
ition) in the arbitrary system to tlie appropriate coordinates in the equinoctial 
system. 
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The transformation from the arbitrary system to the equinoctial system is ex- 
pressed as 


r. « T 


(A-3) 


where 


w page is 
^ quality 


T * R^c-n) R^c.'i R^ca) 


(A-4) 


and r. designates the position vector referred to the arbitrary reference system, 
r^ designates the same position vector referred to either the direct or retrograde 
equinoctial system (depending on the value of the retrograde factor), and where 




cosd imQ 0 

Vmd cjosQ 0 

0 0 1 


(A-5) 


and 




1 0 0 
0 cot>0 sinG 

0 - s'm© C0&© 


(A-6) 


are the matrix representations of the rotation through the angle 6 about the z and 
X axes, respectively. 
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(A-8c) 
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1. V p* ^ C^' 


(A-8d) 


and the transformation matrix is expressed in the equinoctial elements p and q 
as 


T * 




i-pVf 

ipn 

-apl 
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Apcj^X 



(A-9) 

■t 

i 

Ap 

-H 





The rows of this transformation matrix are the components (direction cosines) of 
the T, g, w vectors, respectively, in the arbitrary reference system, i.e.. 
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f - 




iP“l 

-^pl 


(A- 10a) 


The definition of the elements p and q must, of coui'se, be consistent with the 
value of the retrograde factor. 


A-7 


'i ' 
♦ 

>1 


■7 





This is easily demonstrated since 
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A. 3 TIIANSFOHMATION FHOM EC»H lNOCTIAL ELFMFNTS TO POSITION 
AND VELOCITY 

The key to this transfoi'mation is the transformation from oquinot'tial elements 
to position and velocity in tlie equinoctial reference system. The position and 
velocity in any right-hand orthogonal reference system is then obtained by 
inverting the transformation matrix given in Equation (A-7), i.e., 



(A-11) 



(A-12) 


The transformation from equinoctial elements to tlie position and velocity in the 
equinoctial reference system makes use of the mean, eccentric, and true longi- 
tudes, respectively, which are defineil by 


A » X u) IQ 


(A-13i 


F » u. o in. 


(A- 14) 


L » ? -f u) in 


(A- 15) 


where X, u, and f are the mean, eccentric, anil true imomalies. 
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The position and velocity vectors can he expressed as 


Xf ^ (A-16) 

and 


1 *. • C ^ 

Yg (A-17) 


since there is no motion out of the orbital plane. 

Expressions for the coordinates of the position (X, Y) in terms of the tnae long- 
itude follow directly from analytical geometry and are given by 


where 


X = r C0& L 

(A-IS) 

Y = rsiv\l- 

<A-19) 


(A-20) 

1 4- k cos L V\ iin L 



The coordinates of the velocity vector are easily obtained by differentiating the 
expressions for the position coordinates and substituting the following two-body 
relation into the result: 


L = 




(A-21) 


The final results are 


-notCVi » s\r\\S) 


A-22) 


and 


• n 0. ( k cos L) 

Y » - ■ - - - ■ (A-23) 

Vl-h^-k^ 


The position coordinates can be expressed in terms of the eccentric longitude, 
F , using the two-body relations 


r cos(l-<^) « a cos(F- 0) - 


(A-24) 
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0 - cj + Ifl 


(A-26a) 


(A-2(ib) 


The final results are 


X * O, [(l - h*/3) COSF -t- hk/asiriF-kj (A-27) 

and 

Y * tt [(i-k^/9) sirtF hkyflcosF - ] tA- 2 S) 


The velocity coordinates follow by differentiating liquations (A-27) and (A-2J') ami 
substituting the two-bodj' relation 


F a 


0.x no. 

r * r 


(A-2i)i 


yielding 


no. 

r 


[hk^CosF - ( I - h^/8) sim F ] 
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Y = [ ( I - C06 F - hUpVinp] (A-31) 


where the radial distance is expressed as 

r = 0, ( 1 - k cos F - h siv\ F ') 

Equations (A-30) and (A-31) can also be obtained by combining Equations (A-18) 
and (A-19) with Equations (A-27) and (A-28) to yield expressions for cos L and 
sinL. These expressions are substituted into Equations (A-22) and (A-23) to 
yield the final result. 
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A. 4 TRANSFORMATION FROM POSITION AND VELOCITY TO EQUINOCTIAL 
ELEMENTS 

The traiisformation from position and velocity in an arbitrary reference system 
to the equinoctial elements could be obtained by inverting the proper equations 
in Section A. 3. However, appealing directly to the classical two-body problem 
permits a more concise derivation. The semimajor axis is immediately obtained 
by inverting the well known energy integral for the two-body problem which yields 


0. 


1- 

Ifl M. j 


(A-32) 


where “f is the position vector of the satellite in the (x,y,'z) reference system. 
The eccentricity vector is given by 


r 



(f 


(A-33) 


and the unit vector normal to the orbital plane is the normalized angular momen- 
tum vector given by 
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W 




(A-34) 
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In view of Equation (A- 10c) relating the elements p and q to the vector w, it 
follows that 


and 


? 




(A-35) 

ORIGINAE PAGE IS 
OF POOR QUALITY 


^ ’ l+Wjl 


(A-36) 


The elements p and q determined from Equations (A-35) and (A-36) are con- 
sistent with the value of the retrograde factor I . 

The unit vectors f and g may now be computed using Equations (A-lOa) and 
(A-lOb). The equinoctial orbital elements h and k are computed using the 
formulas 


and 


, A 

n a e. • g 


k = e • ^ 


(A-37) 


(A-38) 


The elements h and k are consistent with the vectors f and g with regard to 
the direct and x'etrograde definitions. 
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The remaining element to be computed is the mean longitude, X . First, tlie 
position coordinates X and Y of the satellite relative to the orbital frame g, 
and w are computed from the expressions 


A 

X - r coftU a "r • ^ 


(A-39) 


Y = r sm L « "r • ^ 


(A-40) 


Inverting Equations (A-27) and (A-28\ yields the expi'essions 


COS F a 


k 4- 


- hk/3Y 


(A-41) 


- hk/3X 

s\nF = h +■ n : — 

0. Vl-k'-k*“ 


(A-42) 


which, when substituted into Kepler's equation 


X a F - k s’m F 4 k COS F 


(A-43) 


yields the desired result. 
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A. 5 POISSON BEIACKETS 

In the present application, the Poisson Brackets must be given in terms of the 
equinoctial elements. The results are obtained by direct substitution into the 
previously obtained results of Broucke and Cefola (Reference 33) and are listed 
in Table A-1. 
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Table A-1. Poisson Brackets of Equinoctial Elements 


II 

cf 

- a. 0.^1 

Chik'i » 

"^1^3 

(Ao,h> » 

-h%i| 


- k ps^ 

(Afl, k) * 

-kS4 

(h,(i) = 

- k<^% 

(^o,P) = 

-pSff 

- 

Vipsg 



Ck.cj^ = 





-Ci/aV 


Auxiliary Variables; 


= l/rvCL^ 

Si - a’’ 

“ s^Sj/a + sO 

Sy = s^Sj^/Casa) 


'These expressions are valid for both the direct and retrograde element sets. 
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A.6 PARTIAL DERIVATIVES OF THE EQUINOCTIAL ELEMENTS WITH 
RESPECT TO VELOCITY 


The partial derivatives ba/ir , Sp/a"? , and aq/d"? are obtained directly as 
functions of the equinoctial elements by using the results of Broucke and Cefola 
(Reference 33). However, the expressions for ah/br , ak/br , and b^Q/br 
in terms of the classical orbital elements are not as easily translated into the 
equinoctial elements. To compute these quantities, the following relationship 
(obtained by Broucke (Reference 30)) is used: 



^ ^ ( ^o( > 




bT 

d(Lfi 


(A-44) 


which requires the Poisson Brackets from Table A-1 and the partial derivatives 
of the position vectoi'. To obtain bir/ah and b'r/bk , the following partial der- 
ivatives of X and Y are needed: 


dX 
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kySX 
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(A-45a) 


dX 

dk 


— 4 - (x Y . G) 


(A-45b) 


by 


k/6Y 
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■|(x Y .G) 


(A-45c) 
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bk 
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(A-45d) 
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With these results, the position partial derivatives can be specified, as shown in 
Table A-2. Substitution of the results of Tables A-1 and A-2 into Equation (A-44) 
gives the desired results, which are listed in Table A-3. 
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Table A-3. Partial Derivatives of the Equinoctial 
Elements With Respect to Velocity 
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